A179  741  A  METHOD  TO  OBTAIN  TIHE  DEPENDENT  ELECTRON  EIGENSTATE  1/2 
POPULATIONS  HITH  EL  <U>  MISSION  RESEARCH  CORP  SANTA 
BARBARA  CA  D  H  SOULE  81  JUL  86  MRC-R-1886 
UNCLASSIFIED  DNA-TR-86-282  DNA881-85-C-882S  F/G  28/8  NL 


«.  s 


AD-A179  741 


DNA-TR-86-283 


A  METHOD  TO  OBTAIN  TIME  DEPENDENT  ELECTRON 
EIGENSTATE  POPULATIONS  WITH  ELECTRON  COLLISIONS 
AND  AN  ARBITRARY  RADIATION  FIELD 


D.  H.  Sowle 

Mission  Research  Corporation 

P.  0.  Drawer  719 

Santa  Barbara,  CA  93102*0719 


1  July  1986 


Technical  Report 


CONTRACT  No.  DNA  001-85-C-0035 


Approved  for  public  release; 
distribution  is  unlimited. 


THIS  WORK  WAS  SPONSORED  BY  THE  DEFENSE  NUCLEAR  AGENCY 
UNDER  RDT&E  RMSS  CODE  B322085760  Z99QMXZA00002  H2590D 


Prepared  for 
Director 

DEFENSE  NUCLEAR  AGENCY 
Washington,  DC  20305*1000 


f'DTlC 

fc^EL.£CT 

k  M>R3  0»87 

**  E 


I  . ' .  1, 

r*->js 


k  V1  d 


m 


'.■1 

m 


ft;* 


Destroy  this  report  when  it  is  no  longer  needed.  Do  not  return 
to  sender. 


PLEASE  NOTIFY  THE  DEFENSE  NUCLEAR  AGENCY 
ATTN:  TITL,  WASHINGTON,  DC  20305  1000,  IF  YOUR 
ADDRESS  IS  INCORRECT,  IF  YOU  WISH  IT  DELETED 
FROM  THE  DISTRIBUTION  LIST,  OR  IF  THE  ADDRESSEE 
IS  NO  LONGER  EMPLOYED  BY  YOUR  ORGANIZATION. 


=0R*  SSC.P'V  C-ASS.FCA'CN 

UNCLASSIFIED 


REPORT  DOCUMENTATION  PAGE 


'  3.  RE3TR:C~vE  MAR.C.NGS 


33  DEC_A5S.FCATC.N  DOWNGRADING  SCnEDUL 

N/A  since  Unclassified 


A  RF.RFORMING  ORGANIZATION  REPORT  NUMBER(S) 

MRC-R-1006 


i*.  NAME  OF  PERFORMING  ORGANIZATION 


5.  MONITORING  ORGANiZA*  DN  REPORT  NU.MBERlS) 

DNA-TR-86-283 


6o.  OFFICE  SYMBOL  |  ?i.  NAME  OF  MONITORING  ORGANIZATION 


Mission  Research  Corporation 


:c  AOORESS  Cry,  Sun.  ina  ZIPC odt) 

P.0.  Drawer  719 

Santa  Barbara,  CA  93102-0719 


3a.  NAME  OF  :ljNOiNG  .  SPONSORING 
ORGANIZATION 


3c.  AOORESS  (Cry.  Srar*.  ina  ZIP  Coat) 


(If  wlkioltl 

MRC 


Di  rector 

Defense  Nuclear  Agency 


7b.  AOORESS  (Cry,  Star*,  ina  ZIPC.at) 

Washington,  DC  20305-1000 


BB  OFFICE  SYMBOL  3  PROC-REMENT  NSTRUMENT  OENTIFICATION  NUM8ER 

(if  ipahcioitl 

RAAE/Schwartz  DNA  001-85-C-0035 


'0  SOURCE  OF  =U.NDiNG  NUMBERS 


PROGRAM 


S2715H 


RROjEC* 

NO 

'JO 

Z99QMXZ 

A 

.YORK  _N.' 
ACCES3.CN  NO 

DH  200761 


nauat  Sttunty  Ciisswcitioni 

A  METHOD  TO  OBTAIN  TIME  DEPENDENT  ELECTRON  EIGENSTATE  POPULATIONS  WITH  ELECTRON 
COLLISIONS  AND  AM  ARBITRARY  RADIATION  FIELD 


O  -5RSCNAL  Aw'UOSlSl 

Sowle,  Dave  H. 


13  OAT*  OF  REPORT  Year  Vfonrn.  Oaw  15  RAGE  30 -NT 

860701  118 


'3a  *Y«e  3F  JCRQRT 

Technical 


■5  SUPPLEMENTARY  NOTATION 

This  work  was  sponsored  by  the  Defense  Nuclear  Agency  under  RDT5E  RMS3  Coae  3322035750 
Z99QMXZA00002  H2590D.  _ 


‘3  VwSjcC  _=^'VlSv Confinu#  on  .'ev*rre  r  •'eceiiarv  ma  o enrrrv  ov  oiock  -umoer) 

■r>  Radiative/Colli  si  onal  Transitions  Opaci  ty 
'  Eigenstate  Population  Non-LTE  Opacity 


;9  A3STRACT  Continue  on  r«v«rr«  if  n«c«su/y  and  identify  by  bio <x  numoer) 

\ 

A. computer  module  has  been  written  which  solves  the  time  dependent  electron 
eigenstate  population  problem  under  the  influence  of  electron  collisions  as 
well  as  spontaneous  and  induced  radiative  transitions  to  bound  or  free  states. 
Degree  of  ionization  is  automatically  included.  Electron  velocity  distribution 
must  be  represented  by  a  Maxwellian.  Radiation  is  represented  as  photon  groups, 
called  “features ,  with  flat  central  regions  and  Lorentz-like  tails.  Photon 
features  can  be  superimposed,  so  that  a  rather  general  line  and  continuum  combi¬ 
nation  can  be  represented.  Check  out  examples  discussed  in  the  final  section 
indicate  that  photoionization  may  be  a  more  .important  phenomenon  at  high  altitude 
than  is  generally  appreciated.  v  •  , 
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SECTION  1 
SCOPE 


Many  military  systems,  particularly  SDI  systems,  have  trans- 
nuclear  missions  which  require  successful  propagation  of  electromagnetic 
signals  through  the  atmosphere.  Signals  of  interest  may  be  radio,  radar, 
infra-red  emissions,  optical,  UV,  X-ray,  or  even  atomic  atoms  or  ions. 

Many  of  the  above  wavelength  ranges  are  sensitive  to  the  presence  and 
especially  to  the  structure  of  plasma  or  atomic  and  molecular  emitters  and 
absorbers  produced  or  modified  by  atmospheric  nuclear  explosions.  The 
range  of  yield  and  altitude  of  interest  spans  kilotons  to  tens  of  megatons 
and  ground  surface  to  hundreds  of  kilometers. 

The  atmospheric  test  data  base  is  sparse  except  near  the  sur¬ 
face,  and  lacks  measurements  of  quantities  critical  to  modern  systems  even 
there.  Many  observables  of  interest  cannot  De  obtained  in  undergrouna 
tests.  Consequently  it  is  necessary  to  rely  heavily  on  theoretical  pre¬ 
dictions. 


For  these  purposes  theory  is  required  to  provide  two  distinct 
services.  It  is  required  to  serve  as  a  mechanism  for  interpolating 
between  field  test  events  and  extrapolating  from  them  for  observables 
which  have  been  measured,  and  it  is  required  to  predict  ooservables  which 
have  not  been  measured. 

To  maximize  confidence  in  theoretical  results  for  cases  remote 
from  the  data  base  or  for  observables  not  avail aole  in  the  data  base  it  is 
not  sufficient  to  carry  out  the  most  careful  and  complete  treatment  feas¬ 
ible  on  each  individual  test  event,  or  even  on  several  classes  of  test 


event.  It  is  necessary  to  Go  the  aDove  with  a  single  unified  theoretical 
model,  to  provide  tne  best  assurance  against  an  overlooked  effect  or  a 
misunderstood  transition  in  dominant  physics. 


For  this  reason  an  attempt  to  construct  a  calculation  of  nuclear 
fireball  development  applicable  over  the  entire  range  of  yield  and  alti¬ 
tude  is  underway.  In  addition  to  its  direct  utility  in  predictions  of 
fireball  conditions  for  events  far  from  test  data,  the  major  utility  of 
such  a  model  is  to  provide  predictions  of  the  environment  in  which  early 
time  structure  of  plasma  and  emitters  forms.  This  structure  is  believed 
to  develop  at  later  times  into  the  most  critical  feature  for  many  modern 
systems . 


Nuclear  fireballs  develop  by  processes  of  transport  of  radiant 
and  material  energy.  Basic  to  the  transport  is  the  distribution  of  ion 
species  and  the  electron  eigenstate  populations  within  each  important 
atomic  ion,  for  these  control  ooth  the  chemical  energy  stored  ana  the 
material  opacity.  The  opacity  in  turn  controls  radiation  transport  of 
energy,  under  most  conditions  a  most  important  mechanism.  Therefore  a 
requirement  on  a  unified  theoretical  model  is  the  ability  to  calculate 
atomic  ion  species  distribution  and  electron  eigenstate  distribution 
within  each  atomic  ion  important  either  to  observables  affecting  system 
performance  or  to  radiative  transport  of  energy. 


Depending  on  burst  conditions  ano  the  system  requirements,  the 
calculation  must  be  reliable  for  times  from  microseconds  to  minutes,  for 
temperatures  from  kilovolts  to  a  few  hundred  degrees  Kelvin,  and  for  elec¬ 
tron  densities  from,  say,  106/cm3  or  less  to  10  22/cm3  or  more.  Further, 
it  must  handle  radiation  fields  from  sparse  line  spectra  to  Planckian.  In 
particular,  it  must  handle  conditions  which  span  the  range  from  extreme 
non- thermodynamic  equilibrium  to  local  thermodynamic  equilibrium. 
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This  report  describes  a  method  developed  to  satisfy  this  need 
tes  monatomic  ion  SDecies  populations  and  electron  eigenstate 


SECTION  2 
REQUIREMENTS 


2.1  SYSTEMS. 

At  a  level  this  deep  in  theoretical  modelling  the  only  system 
requirement  laid  on  the  method  is  that  it  be  capable  of  handling  chemical 
elements  important  to  systems  (emitters,  absorbers,  etc.)  but  not  impor¬ 
tant  to  energy  or  material  transport.  This  is  simply  accomplished,  by 
allowing  for  mixtures  of  a  large  number  of  arbitrary  elements  to  be  pres¬ 
ent  at  any  point  in  space.  The  price  paid  for  this  capability  is  some 
storage  on  nonvirtual  memory  computers  and,  presumably,  almost  nil  on 
virtual  menory  machines. 

2.2  PHYSICS. 

Consider  a  chemical  element  immersed  in  an  field  consisting  of 
hot  ions,  electrons,  and  photons.  It  is  necessary  to  find  the  rate  at 
which  energy  is  being  absorbed  or  emitted  by  this  element.  To  do  so  one 
must  know  the  distribution  of  ionization,  or  charge,  states  of  the  ele¬ 
ment,  the  distribution  of  electrons  in  the  eigenstates  of  each  charge 
state,  and  the  rate  of  collisional  and  radiative  processes  from  each 
eigenstate  to  all  others.  Perhaps  the  simplest  way  to  think  of  the  prob¬ 
lem  is  to  consider  a  single  array  of  eigenstates  for  an  element;  that  is, 
consider  the  ground  state  of  the  lowest  charge  state  (neutral  or  negative 
ion)  as  the  base  state  and  all  other  states,  including  higher  ionization 
states,  as  excited  states.  In  fact,  this  has  turned  out  to  be  the  best 
way  to  accomplish  the  numerical  solution.  In  the  case  of  hydrogen  one  has 


the  ground  state  of  the  neutral  atom,  call  it  H°(l),  and  its  excited 
states  plus  a  single  eigenstate  corresponding  to  a  bare  proton,  H^l).  In 
the  case  of  oxygen  one  has,  at  present,  the  ground  state  of  an  oxygen  atom 
with  attached  electron,  01-(I),  followed  by  the  ground  state  neutral, 

0*(1)  at  about  an  eV  excitation  energy,  then  all  the  excited  states  of  0°, 
followed  by  all  the  states  of  01,  through  08(1),  the  bare  oxygen  nucleus. 

The  basic  problem  then  is  to  find  the  electron  distribution 
among  eigenstates  of  an  element.  Once  that  is  done,  charge  distribution, 
energetics,  and  opacity  all  follow  easily. 

The  processes  which  act  to  redistribute  electrons  among  the 
eigenstates  include  electron  collisions,  radiatively  induced  transitions, 
spontaneous  radiation  and  heavy  ion  collisions.  Heavy  ion  collisions  are 
less  important  than  the  others  and  have  been  neglected,  although  the 
numerical  solution  is  so  organized  that  heavy  ion  processes  can  be  intro¬ 
duced  in  a  very  straightforward  manner  if  they  ever  appear  to  be  neces¬ 
sary.  Another  restriction,  which 'presently  appears  to  be  more  important, 
is  that  di-electronic  recombination  and  its  inverse,  Auger  type  processes, 
are  neglected.  Di-electronic  recombination  sometimes  is  the  most  impor¬ 
tant  --ecomDi nation  Miechanism  under  conditions  of  high  electron  temperature 
and  low  electron  density.  But  under  those  conditions  very  little  recombi¬ 
nations  occurs;  this  is  the  justification  for  reducing  the  complexity  of 
the  initial  code  version  by  neglecting  such  processes. 

The  processes  treated  are:  bound-bound  electron  collisions 
(upward  and  downward  in  excitation  energy) ,  bound-free  electron  collisions 
along  with  the  inverse  three  body  recombination,  bound-bound  radiation 
induced  transitions  up  and  down,  bound-free  radiation  induced  transitions 
'ionization  and  recombination),  and  bound-bound  spontaneous  emission  (down 
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For  most  rates  a  factor  of  two  to  three  is  the  limit  of  current 
knowledge;  for  some  important  ones  the  situation  is  better  but  for  many  it 
is  even  worse.  A  factor  of  three  error  in  a  single  rate  coefficient  can 
result  in  as  much  as  a  factor  of  three  error  in  population,  while  an 
unfortunate  coincidence  involving  such  errors  in  several  rate  coefficients 
can  result  in  even  greater  error  in  population. 

Uncertainty  in  the  correct  value  of  rate  coefficients  is  a 
fundamental  limitation  on  accuracy  of  theory.  In  a  situation  as  complex 
as  this,  the  only  satisfactory  method  of  approach  is  sensitivity  analysis, 
where  the  effects  of  modification  in  key  rate  coefficients  are  assessed  by 
inserting  modified  values  into  the  code  and  comparing  results. 


2.3 


MATHEMATICS . 


The  major  mathematical  requirement  on  the  method  is  generality, 
due  to  the  extremely  broad  range  of  conditions  over  which  the  method  is  to 
be  applied.  This  means  that  one  is  not  allowed  to  use  most  normal  methods 
designed  for  special  limiting  cases;  LTE,  coronal,  etc.,  but  must  adopt  a 
general  method  which  is  capable  of  handling  both  those  limiting  cases  and 
intermediate  cases.  In  short,  the  equations  must  be  solved  directly. 


Given  that  important  physical  parameters,  especially  rate  coef¬ 
ficients  but  also  energy  levels  in  some  cases,  have  substantial  uncertain¬ 
ties,  the  only  firm  mathematical  accuracy  requirement  on  the  method  is 
that  it  yield  answers  accurate  to  better  than,  say,  a  factor  of  three  or 
so.  Unfortunately  it  is  difficult  if  not  impossible  to  analyze  the  final 
effect  of  a  mathematical  approximation  buried  deep  in  such  a  complex  set 
of  equations  as  those  which  must  be  solved,  so  it  is  most  unwise  to  allow 
any  unnecessary  inaccuracies  in  the  math.  Certainly  it  is  almost  always 
easier  to  do  the  math  correctly  than  to  analyze  the  implication  of  an 
approximation.  Nevertheless  it  has  been  necessary  to  give  up  mathematical 


rigor  at  a  few  points  in  the  method  to  meet  finite  computer  time  require¬ 
ments.  The  impact  of  such  approximations  on  the  final  results  have  been 
assessed  by  comparison  of  these  calculations  with  results  obtained  from 
more  precise  methods.  Enough  precise  results  are  available  to  provide 
reasonable  assurance  of  the  mathematical  soundness  of  our  method. 

2.4  COMPUTER. 

Computer  requirements  are  that  the  coded  implementation  of  the 
method  be  transportable  among  mainframes,  fit  within  memory  available  on  a 
Cray-class  computer  and  that  this  portion  of  the  code  require  at  most  a 
few  hours  of  CP  time  on  such  equipment  to  complete  calculation  of  a  fire¬ 
ball  history. 

The  transportability  requirement  creates  some  tedious  but 
straightforward  limitations  on  code  and  data  structure,  mainly  due  to  the 
fact  that  Crays  do  not  have  virtual  memory  operating  systems. 

The  memory  requirement  does  not  appear  to  be  limiting. 

However,  the  CP  time  requirement  has  been  the  major  driving 
force  in  almost  every  aspect  of  design  and  implementation  of  the  method. 
The  time  allowed  for  any  method  implemented  to  solve  the  eigenstate  popu¬ 
lation  equations  can  be  estimated  as  follows.  As  an  average,  take  M 
eigenstates  per  ionization  state,  Z  ionization  states  simultaneously  pres¬ 
ent  per  element,  E  elements  present  per  spacial  zone,  C  space  zones,  and  T 
time  increments  per  problem.  Then  the  number  of  eigenstate  solutions 
necessary  to  complete  a  problem  is 

N  =  MZECT  .  (1) 

A  reasonable  minimum  number  of  eigenstates  necessary  to  get  a 
decent  value  of  opacity  is,  say, 
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It  is  difficult  to  imagine  a  scheme  so  wise  that  it  can  solve  a 
problem  without  computing  for  3  or  4  ionization  states  simultaneously 
present  under  severe  conditions,  but  perhaps  only  half  the  space  cells 
will  have  severe  conditions,  take 

Z  =  2 

The  number  of  elements  per  cell  is  at  least  2  (N  and  0)  and  more 
1 ikely  to  average  3, 

E  =  3. 

Space  zones  depend  on  the  dimensionality  of  the  problem,  let  D 
be  the  number  of  space  dimensions,  then  C  is  at  least  100  per  dimension. 


Finally,  one  expects  a  problem  to  be  completed  in  about  a  thou¬ 
sand  time  steps,  so  take 

T  =  1000. 

Then 

N  =  6  x  104+2° 

and  if  the  total  time  spent  on  the  method  is  to  be  H  hours,  the  time 
allowed  per  eigenstate  solution  is 


t  =  1.7  x  10- 


5  3600  H  _  6  x  10- 2  H 
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Given  that  H  is  to  be  of  order  1,  we  find  that  for  a  three  dimensional 
code 

t  =  6  x  10" 8  sec 

This  is  hopeless  for  foreseeable  equipment.  To  have  a  reason¬ 
able  chance,  we  must  limit  ourselves  to  one  space  dimension  so  that 

t  =  6  x  lO"4  sec  =  0.6  msec 

This  translates  into  12  msec  to  solve  all  eigenstates  and  charge 
states  for  a  given  element  in  one  spacial  zone  per  time  cycle.  Not  much 
time.  To  date  we  are  within  a  factor  of  two  to  ten  of  the  goal,  depending 
an  the  complexity  of  the  radiation  field,  *or  the  MRC  ELXSI  computer  va  4 
to  5  MIP  machine).  All  things  being  as  advertised,  this  should  place  us 
at  or  better  than  the  goal  on  a  CRAY. 


METHOD  ADOPTED  FOR  EIGENSTATE  SOLUTION 


3.1  SUMMARY  OF  METHOD. 

(a)  Transition  rates  for  all  processes  are  calculated  connecting 
states  of  Z  to  al 1  eigenstates  of  Z-l,  Z,  and  Z+l. 

(t>)  Eigenstates  for  all  ion  charge  states  of  a  given  element  are 
arranged  in  a  linear  array  ordered  on  increasing  state 
energy  referenced  to  the  lowest  energy  state  of  the  element 
(ground  state  neutral  or  ground  state  negative  ion).  Call 
the  total  number  of  equations  in  this  array  L. 

(c)  The  total  state  array  is  sorted  into  two  linear  arrays,  one 
set  to  be  solved  with  a  quasi-steady  approximation  and  the 
other  to  be  solved  as  differential  equations. 

(d)  These  two  sets  of  equations  are  each  further  broken  down 
into  sets  of  closely  coupled  equations,  for  a  total  of  M 
sets,  1  <  M<  L. 

(e)  Each  of  the  M  sets  is  solved  by  Gauss  elimination  using 
initial  conditions  derived  from  the  original  linear  array  of 


(f)  Sets  initially  containing  a  major  fraction  of  the  total  ele¬ 
mental  electron  population  are  renormalized  to  their  initial 
total  populations.  Sets  initially  containing  a  minor  frac¬ 
tion  of  the  total  elemental  electron  population  are  allowed 
major  population  changes,  but  minor  changes  are  strongly 
damped. 

(g)  Effective  set-to-set  transition  rates  are  calculated  based 
on  individual  eigenstate  rates  and  populations. 

(h)  The  pair  of  sets  with  the  fastest  interset  rate  is  solved  as 
a  pair  of  coupled  differential  equations  for  total  set  popu¬ 
lation,  with  relative  eigenstate  population  unchanged. 

(i)  The  resultant  sets  from  step  (h)  are  amalgamatea  into  a 
single  set,  reducing  the  number  of  sets  by  1. 

(j)  Steps  (h)  through  (i)  are  repeated  until  only  one  set, 
encompassing  all  eigenstates  of  all  ionization  states, 
remains . 

(k)  Because  absolute  normalization  was  lost  in  step  (f),  the 
final  set  is  renormalized  to  conserve  atomic  ions. 

3.2  SUMMARY  OF  METHOD  RATIONALE. 

The  overriding  objective  in  developing  the  method  was  to  mini¬ 
mize  CP  time  without  too  much  loss  in  accuracy.  The  principal  tactic  has 
been  to  maximize  allowable  time  step.  The  process  consisted  of  recogni¬ 
tion  of  a  problem,  curing  it,  recognition  of  problems  created  by  the  cure, 
curing  them  etc.  until  no  new  problems  appeared.  The  result  is  cumbersome 
logically,  but  it  is  quite  rugged,  adequately  accurate,  and  relatively 
fast  in  execution. 


The  first  step  is  to  decide  how  many  ionization  states  and 
eigenstates  per  ion  must  be  calculated.  This  is  done  before  any  rates  are 
calculated.  Rate  calculation  is  expensive;  the  order  of  102  rates  and 
inverses  must  be  calculated  for  each  eigenstate. 

All  rates  connecting  each  eigenstate  of  ion  Z  to  all  other 
eigenstates  present  of  ions  Z-l,  Z,  and  Z+l  are  then  calculated. 

To  make  it  automatic  that  eigenstate  equations  in  adjacent  ioni¬ 
zation  states  be  grouped  together  for  solution  where  appropriate,  eigen¬ 
state  data  for  all  values  of  Z  present  are  stacked  in  a  single  linear 
array  indexed  on  increasing  eigenenergy.  This  arrangement  turns  out  to 
save  some  storage  and  simplify  some  logic,  as  added  bonuses. 

The  method  is  required  to  work  for  situations  where  some  or  all 
states  approach  LTE.  This  guarantees  a  severe  stiff  differential  equation 
problem.  Stiff  equations  can  be  solved  assuming  them  to  be  quasi-steady. 

A  criterion  to  distinguish  between  stiff  and  "limp"  behavior  is  straight¬ 
forward;  but  some  efficiency  must  be  sacrificed  here.  We  must  allow  equa¬ 
tions  to  go  from  stiff  to  "limp"  as  well  as  from  limp  to  stiff  and  unfor¬ 
tunately  criteria  which  identify  every  stiff  equation,  once  met,  force 
equations  to  remain  stiff  forever.  Accordingly  a  less  efficient  criterion 
is  applied  and  the  linear  array  of  eigenstate  equations  is  sorted  into  two 
arrays  or  sets,  the  QS  set  to  be  solved  under  a  quasi-steady  assumption 
and  the  DE  set  to  be  solved  as  ordinary  differential  equations. 

Under  circumstances  of  interest,  normalization  problems  exist  if 
the  two  sets  are  solved  in  a  straightforward  manner.  Some  or  all  of  these 
problems  were  created  by  breaking  the  equations  into  two  sets.  To  address 
these  problems,  each  of  the  two  sets,  QS  and  DE,  is  further  broken  down 
into  sets  of  closely  coupled  equations. 


Each  QS  set  is  then  solved  by  Gauss  elimination  separately, 
using  initial  conditions  from  all  other  sets  as  input,  and  each  DE  set  is 
similarly  solved  separately  by  a  different  version  of  Gauss  elimination. 

The  resultant  set  electron  populations  are  incompatible  in  gen¬ 
eral,  since  changes  in  other  set  populations  were  neglected  when  finding 
eigenstate  population  in  a  given  set.  The  basic  way  adopted  to  handle 
this  is  to  renormalize  each  set  back  to  its  initial  population,  so  that 
only  the  electron  distribution  within  each  set  has  been  changed,  calculate 
effective  rates  between  the  sets,  then  solve  differential  equations  for 
set  populations. 

The  situation  at  this  point  is  similar  to  that  at  the  beginning 
of  the  process,  a  number  of  linear  differential  equations  for  set  (rather 
than  eigenstate)  populations  must  be  solved.  If  we  were  guaranteed  that 
the  number  of  sets  were  less  than  the  number  of  initial  equations  we  could 
simply  repeat  the  process  (sort  into  QS  and  DE,  sort  into  closely  coupled 
equations,  solve)  iteratively  until  only  one  set,  including  all  eigen¬ 
states,  remained.  Unfortunately,  it  is  not  true  that  under  all  circum¬ 
stances  each  such  iteration  would  reduce  the  number  of  sets. 

The  technique  actually  used  is  to  analytically  solve  only  for 
populations  of  the  two  most  closely  coupled  sets,  neglecting  effects  of 
the  others.  These  two  sets  are  amalgamated  into  a  single  set,  rates  among 
the  reduced  number  of  sets  are  adjusted  and  the  process  repeated  until 
only  one  set  remains. 

This  time  splitting  technique  for  set  populations  has  its  own 
version  of  the  stiff  equation  problem.  It  is  always  far  less  severe  than 
that  of  the  original  eigenstate  set  and  appears  only  to  affect  sets  with 
minor  fractions  of  total  electron  population,  but  it  is  bao  enough  to 
require  fixing.  The  fix  adopted  is  empirical  rather  than  rigorous  and  is 


targeted  against  the  specific  problem  encountered.  Because  the  time 
splitting  scheme  locks  set  populations  together  by  pairs,  a  low  population 
set  strongly  influenced  by  two  high  population  sets  and  requiring  a  major 
population  adjustment  can  become  locked  into  a  position  where  the  adjust¬ 
ment  is  impossible.  To  avoid  this  eventuality,  at  the  end  of  the  Gauss 
elimination  procedure,  the  Gauss  answer  for  set  population  is  accepted  for 
minor  population  sets  with  major  population  changes,  but  all  sets  with 
major  population  are  renormalized  as  are  minor  population  sets  with  minor 
population  changes. 

With  this  modification  the  process  of  sorting  equations  into 
sets,  solving  for  relative  populations  within  the  sets  by  Gauss  elimina¬ 
tion,  then  solving  for  set  population  by  time  splitting  works  adequately. 
It  is  only  necessary  to  add  a  final,  overall,  normalization  since  we  have 
allowed  the  absolute  normal ization  to  be  broken  occasionally  after  the 
Gauss  elimination  step. 

3.3  MAJOR  ELEMENTS  OF  THE  METHOD. 

3.3.1  Adding  Ionization  States. 

A  time  cycle  begins  for  some  chemical  element  with  one  or  more 
contiguous  ionization  states,  Z,  present,  and  with  a  distribution  of  elec¬ 
trons  in  their  eigenstates.  Circumstances  may  require  addition  of  another 
state  of  ionization  one  unit  charge  above  the  highest  present  or  one  unit 
below  the  lowest  present.  In  preparation  for  this,  an  ion  state  above  the 
highest  and  one  below  the  lowest  present  are  added  without  populations,  if 
such  ion  states  exist  for  the  element.  At  this  point  the  maximum  number 
of  charge  states  which  may  possibly  be  considered  has  been  determined. 


3.3.2 


Deleting  Ionization  States. 


It  is  important  to  avoid  calculating  and  processing  unnecessary 
rates,  to  save  time.  Before  any  rates  are  calculated,  "temperatures"  are 
defined  for  all  processes  which  might  control  creation  or  destruction  of 
charge  states.  At  present  there  are  three  such  temperatures,  one  is  elec¬ 
tron  temperature,  one  is  the  temperature  of  that  Planck  spectrum  which 
contains  the  same  total  radiant  flux  as  the  actual  (line  or  continuum) 
spectrum,  and  one  is  that  of  a  grey  body  spectrum  with  the  same  mean  pho¬ 
ton  energy  as  the  actual  spectrum. 

The  highest  density  ion  species,  D(Zm),  is  then  determinea. 

To  eliminate  the  lowest  state  of  ionization,  Z,  the  equilibrium 
'-atio,  R  ,  of  Z  to  Z+l  is  estimated  using  the  colaest  of  the  three  temper¬ 
atures,  then  it  must  be  true  that  both 

D(Z)  .  r  D(Zm) 


and  either 


D(  Z+l )  <  f  0(Zm) 


or 


R 


e 


<  R 


q 


where  the  current  value  of  r  is  10-3  and  f  is  currently  0.1  if  D ( Z ) =0  or 
0.05  if  D(Z)  is  nonzero  (to  make  it  more  difficult  to  eliminate  an 
existing  ion  state  tnan  one  just  added),  "he  value  of  R  used  is  10*'. 


If  an  ion  state  is  eliminated,  its  density  is  added  to  Z+l 
(without  modifying  eigenstate  distribution  in  Z+l)  and  the  process  re¬ 
peated  until  the  test  (2)  fails. 

To  eliminate  the  highest  state  of  ionization,  Z,  the  equilibrium 
ratio,  R  ,  of  Z  to  Z-l  is  estimated  using  the  hottest  of  the  three 
temperatures  and  the  sequence  of  tests  (2)  applied  to  it.  Again  this  test 
is  repeated  until  it  fails. 

At  this  point  only  those  ion  states  with  significant  popula¬ 
tions,  plus  those  which  may  gain  significant  populations  in  this  cycle  are 
present . 

3.3.3  Adding  an  Artificial  Eigenstate. 

There  are  two  reasons  for  limiting  the  numDer  of  eigenstates  per 
;onizat’on  state  to  aoout  10,  wnicn  implies  principal  quantum  njinoer 
limited  to  3  or  4.  The  first  reason  is  that  reliaole  data  ’s  generally 
very  incomplete  at  higher  quantum  numbers  for  elements  of  interest  so  that 
a  more  complete  representation  would  be  interspersed  lioe^ally  with  e'tner 
gaps  or  guesswork.  The  second  reason  is  that  CP  restrictions  do  not  allow 
a  large  number  of  eigenstates  to  oe  represented.  The  lumpe1”  of  state- te¬ 
state  rate  coefficients  which  must  be  calculated  is  oroportional  to  the 
square  of  the  number  of  states  represented,  and  time  spent  in  the  Gauss 
elimination  routine  is  proportional  to  the  cube  of  the  number  of  states. 

Nevertheless,  there  may  be  tens  to  thousands  of  eig  nstates 
between  the  uppermost  state  represented  and  the  ionization  continuum.  The 
higher  states  must  be  taken  into  account  if  correct  ionization  and  recom¬ 
bination  rates  to  the  interesting  lower  states  are  to  be  obtained.  This 
is  accomplished  by  creating  an  artificial  eigenstate  to  reo^esent  the 
upper  states  not  modeled. 


To  create  such  a  state  one  first  finds  the  number  of  states 
lying  between  the  uppermost  state  modeled  and  the  ionization  continuum. 
This  is  determined  by  the  amount  by  which  the  ionization  potential  is 
depressed  due  to  polarization  effects  of  the  plasma.  From  reference  1, 
pp.  137-140,  the  energy  of  this  depression  is  roughly 

Ed  =  (Z+l)  e2/^  (3) 

where  Z  is  the  degree  of  ionization,  e  (esu)  is  electron  charge,  and  pq 
is  Debye  length, 

PD  =  [kT/(4ne2Ne(l+<Z2>/<Z>))  j1/2  (4) 

with  k  Boltzmann's  constant,  T  the  electron  temperature,  and  electron 
density. 

Assuming  the  upper  levels  to  be  hydrogenic,  the  principal  quan¬ 
tum  number  of  the  level  which  lies  at  the  ionization  continuum  is  given  by 

Eh»*1)2/s2  .  E„  (5) 

or 

sM  ■  (Z+l)  Xfo  . 

where  is  the  ionization  potential  of  hydrogen  in  vacuo.  States  above 
Sm  are  ignored.  If  the  principal  quantum  number  of  the  highest  energy 
state  represented  is  ,  then  the  artificial  state  must  represent  all 
states  from  Sfj  to  S^. 

Parameters  to  represent  the  artificial  state  must  be  calcu¬ 
lated.  The  basic  parameters  are  calculated  as  simple  integrals  weighted 


by  the  degeneracy  of  the  upper  states.  For  hydrogenic  states  the  degener¬ 
acy  of  a  state  with  principal  quantum  number  n  is 


gn  =  2n 2  .  (6) 

Then  the  degeneracy  of  the  artificial  state,  g A,  is  the  sum  of 
degeneracies  of  its  constituents  multiplied  by  the  degeneracy  of  the 
ground  state  of  the  next  higher  state  of  ionization, 

9A  *  9Z+1  /  2"2dn  =  §  (S3-S3)(2S++l)(2l++l)  (7) 

where  S+  is  the  spin  quantum  number  of  the  ground  state  of  Z+l  and  L+  is 
the  orbital  quantum  number  of  Z+l. 

The  mean  energy  of  the  artificial  state,  EA,  is 
E A  =  V  f  (l-Q2/n2)2n2dn/ ;  2n2dn 

*  V(1-3Q2(Sm-Sn)/(Sm3-Sn3))]  (3) 


where  V  is  the  vacunn  ionization  potential  of  the  ground  state  ion  Z  and  Q 
is  the  ground  state  of  Z  principal  quantum  number.  The  factor  Q2  scales 
the  ionization  potential,  V,  to  its  equivalent  hydrogenic  value. 

Similarly,  the  mean  principal  quantum  number  of  the  artificial 
state,  QA,  is  given  by 

Qa  -  /  n2n 2dn/ /  2n2dn  =  1  ( S^-S^ ) / ( S^-S^ )  .  (9) 


Parameters  for  Stark  broadening  and  for  rate  coefficient  calcu¬ 
lations  must  also  be  generated.  The  Stark  broadening  term  for  the  state 
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is  estimated  from  the  lead  term  in  equation  (46)  and  the  crucial  oscilla¬ 
tion  strength  calculations  are  discussed  in  section  5.7.2. 


3.3.4  Deleting  Eigenstates. 


3. 3. 4.1  Before  Calculations  Commence.  In  a  situation  where  the  hottest 


effective  temperature  (see  3.3.2)  is  small  compared  to  the  excitation 
energy  of  some  eigenstate  it  may  be  appropriate  to  save  some  calculation 
time  by  neglecting  that  eigenstate.  Accordingly,  after  unnecessary  ions 
have  been  removed  as  described  in  section  3.3.2  the  three  following  tests 
are  made,  from  highest  excitation  state  for  which  data  exists  to  the 
ground  state  of  the  highest  ionization  state  present. 


Is  the  equilibrium  ratio  of  ion  Z+l  to  ion  Z  greater  than  about 


Is  the  fractional  population  of  eigenstate  s  greater  than 


10' 3? 


Is  the  energy  of  state  s  less  than  ten  times  the  hottest  temper¬ 


ature? 


If  any  of  the  three  tests  is  met,  then  eigenstate  s  and  all 
lower  eigenstates  are  accepted.  If  all  tests  fail,  s  is  deleted  and  the 
tests  are  applied  to  s-1. 


The  principal  virtue  of  this  procedure  is  to  eliminate  upper 
states  from  consideration  under  benign  conditions.  For  example,  undis¬ 
turbed  cells  at  the  beginning  of  a  problem,  or  material  late  in  a  problem 
which  has  recombined  and  fallen  to  the  ground  state. 


3. 3. 4. 2  During  Calculations.  If  the  value  of  Ep  from  (3)  is  so  large 
that  the  energy  from  some  states  s,  Es,  satisfies 
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then  these  states  are  removed. 

There  are  two  physics  problems  associated  with  this  procedure. 
The  minor  one  is  that  all  eigenstates  lying  above  the  effective  ionization 
potential 

A  =  V-Ed  (11) 

still  exist,  even  though  they  lie  above  the  ground  state  of  Z+l.  In 
principle,  account  should  be  taken  of  these  states.  One  straightforward 
procedure  might  be  to  add  the  combined  oscillator  strength  to  these  states 
from  state  s  to  that  from  state  s  to  the  ground  state  of  Z+l,  on  the 
assumption  that  electrons  raised  to  states  above  the  ground  state  of  Z+l 
will  be  lost  (of  course,  they  might  even  lie  above  some  excited  state  of 
Z+l,  to  complicate  matters).  A  procedure  such  as  this  has  not  been  imple¬ 
mented  because  the  numerical  effect  doesn't  seem  important  under  condi¬ 
tions  where  our  treatment  is  credible  and  more  serious  problems  exist 
where  the  treatment  becomes  questionable  due  to  the  second  physics  prob¬ 
lem,  as  fol lows. 

If  electrons  become  so  closely  packed  that  their  Debye  spheres 
begin  to  overlap  (at  1-eV  this  would  be  Ng  >  1020)  the  potential  energy 
wells  of  neighboring  ions  begin  to  overlap;  the  plasma  begins  to  show 
liquid  or  solid  properties!  The  concept  of  an  isolated  ionic  system  of 
eigenstates  begins  to  fail,  electrons  no  longer  have  all  space  available 
for  orbits,  they  must  find  channels  between  ions.  We  have  not  addressed 
this  problem  and  hope  never  to  be  required  to  do  so. 

What  has  been  done  is  entirely  ad  hoc.  The  last  4  eigenstates 
are  retained  independent  of  the  outcome  of  the  three  tests  discussed  in 
the  foregoing. 


3.3.5 


Stiff  Equation  Problem. 


To  find  the  populations  of  al 1  M  eigenstates  of  (all  ionization 
states  of)  some  element,  the  solution  is  required  to  the  set  of  equations: 

M  M 

dNs/dt  =  T$  oksNk  -  N$  +  Bs]  (12) 

for  s  from  1  to  M. 

In  (12),  all  terms  are  inherently  non  negative,  a-jj  is  the 
rate  of  transfer  from  state  i  to  state  j  due  to  all  causes,  per  unit  popu¬ 
lation  of  state  i  and  per  unit  time,  Ni  is  population  of  state  i  and  the 
terms  Ts  and  Bs  refer  to  sources  and  sinks  of  population  outside  the 
set  of  M  equations. 

Equation  sets  like  (12)  are  infamous  for  having  so  called 
"stiff"  equation  difficulties  when  any  straightforward  method  of  solution 
is  applied.  This  occurs  whenever  for  some  value  of  s,  the  input  terms 
tend  to  cancel  the  output  terms,  driving  the  time  step  toward  zero  while 
Ns  is,  in  fact,  nearly  constant  in  time.  In  our  case  such  a  problem  is 
guaranteed  since  we  wish  to  handle  cases  where  equilibrium  is  approached, 
where  the  inputs  and  outputs  effectively  cancel  for  all  s.  Given  that  a 
one  dimensional  radiation  transport  problem  will  typically  consist  of 
thousands  of  radiation  and  hydro  time  steps,  hundreds  of  cells,  and  sev¬ 
eral  elements  per  cell  (N,0,  etc.)  with  perhaps  40  eigenstates  per  element 
(say  10  states  each  for  4  ionization  states),  a  stiff  equation  problem 
would  be  catastrophic. 

To  avoid  stiff  equation  difficulties  the  M  equations  are  sorted 
into  a  set  of  "Quasi-steady"  equations  (hereafter  referred  to  as  "QS")  and 
a  set  of  "Differential"  equations  (hereafter  referred  to  as  "DE")  which 
are  treated  separately.  This  separation  is  currently  based  on  the  condi¬ 
tion: 


where  At  is  the  time  step,  arbitrary  for  purposes  of  this  test,  and  G  is 
some  sizable  number;  G  =  e3  =  20.1  is  the  value  currently  in  use,  but  the 
precise  value  of  G  does  not  appear  to  be  critical. 

This  test  is  made  independently  each  time  the  set  is  to  be 
advanced,  so  a  given  equation  is  able  to  be  advanced  as  QS  or  DE  as  appro¬ 
priate  to  current  conditions;  no  equation  becomes  frozen  into  one  mode  or 

the  other. 

3.3.6  Normalization  Problems. 

Three  classes  of  normalization  problems  arise  when  solving  these 
equations  by  the  chosen  method. 

3. 3.6.1  Inherent  Numerics  Problem.  The  first  class  of  normalization 
problems  is  inherent  in  eq  (12).  To  demonstrate  it,  sum  (12)  over  all  s 
to  obtain 

d(ENs)/dt  =  ITS  -  INSBS  (14) 

where  the  terms  containing  a' s  have  cancelled.  If  the  set  were  isolated 
as  it  would  have  been  had  we  not  separated  the  equations  into  two  classes, 

then  Ts  and  Bs  would  be  zero  for  all  s  and,  ideally,  eNs  would  be 

constant.  Almost  any  method  of  numerical  solution  will  lose  this  preci¬ 
sion  but  that  is  easy  to  correct;  simply  slid  Ns  before  and  after  solu¬ 
tion  then  multiply  the  new  values  by  the  ratio  to  force  conservation  of 
total  number. 


3. 3. 6. 2  Set  Separation  Problem.  The  second  class  of  normal i zation  prob¬ 
lem  is  caused  by  separation  of  the  equations  into  sets  to  be  treated  dif¬ 
ferently.  Now  Ts  and  Bs  are  no  longer  zero  in  general.  It  Decomes 
possible  to  lose  or  gain  population  excessively  from  one  set  to  another. 
The  problem  is  even  deeper  than  it  appears,  for  it  is  not  only  possible 
but  fairly  common  for  a  subset  of  DE  equations  to  be  coupled  more  tightly 
to  some  subset  of  QS  equations  than  it  is  to  the  remaining  DEs,  and  the 
reverse  also  occurs.  So  it  is  not  necessarily  true  that  a  common  normali¬ 
zation  factor,  however  derived,  will  be  appropriate  to  all  DEs  or  all  QSs. 

After  the  separation  into  DE  and  QS  a  further  separation  into 
subsets  is  made  to  locate  equations  which  are  significantly  more  closely 
coupled  to  each  other  than  to  others.  The  precise  method  for  accomplish¬ 
ing  this  is  discussed  in  section  3.3.7. 

Depending  on  circumstances,  and  the  separation  algorithm,  after 
t.nis  process  there  are  a  nuinoer  of  distinct  equation  sets  wnich  may  vary 
from  I  (all  are  DE  or  all  are  QS  ana  all  are  reasonaDly  coupled  together] 
to  M  (none  tightly  coupled  to  another  in  the  same  class).  The  next  steD 
in  the  solution  is  to  collect  the  terms  in  equation  (12)  for  each  subset. 
Ts  ana  8S  thus  are  constructed  from  ii<;S  ana 

Ts  -  "  'V  Vs  •  3s  =  ’sk1  15i 

k*  k’ 

where  k'  does  not  belong  to  the  set  to  which  s  belongs.  After  solution, 
most  sets  are  renormalized  to  a  value  such  that  total  population  within 
the  set  is  close  to  its  original  value  (see  section  3.3.8).  Thus  the  step 
mainly  redistributes  electrons  among  the  eigenstates  of  the  subset,  with¬ 
out  allowing  the  total  population  of  the  subset  to  change. 


Next,  effective  a  s  are  constructed  connecting  each  set  m  to  ail 
other  sets  n 


nm 


J  Ns  l  ask7  J  N. 
s  k  s 


;i6) 


where  s  belongs  to  set  n  and  k  belongs  to  set  m. 

Then  values  of  n‘  and  m‘  are  found  for  which  «  .  ,  +  a.  .  is  a 

n  m  in  n 

maximum.  These  two  sets  are  then  solved  as  a  differential  equation  inde¬ 
pendent  of  all  other  sets. 


N  =  a  N  -  a  N 
n  ran  m  nm  n 


yt»«)  .  yt»E  *  Tyu„^nHHi  a?) 

where 

T  ?  Nn(t)+N;n(  t) 

This  solution  preserves  normal i zation  and  finds  a  new  value  for 
the  relative  populations  of  eigenstates  in  the  sets  n'  and  m‘ . 

Sets  n‘  and  m'  are  then  amalgamated  into  a  single  set,  set-to- 
set  ct*  s  readjusted,  and  the  process  repeated  until  only  one  set,  compris¬ 
ing  all  eigenstates,  remains.  The  solution  is  then  complete. 


i 


This  rather  peculiar  form  of  time  splitting  has  the  advantage  of 
giving  priority  to  the  fastest  reactions  and  of  approaching  the  correct 
steady  state  solution.  However,  if  At  is  too  large,  it  will  approach  the 
steady  state  solution  too  slowly;  its  behavior  is  similar  to  an  implicit 
difference  scheme  in  this  respect. 

3. 3. 6. 3  Gauss  Elimination  Problem.  The  third  and  final  class  of  normal¬ 
ization  problems  arises  from  the  particular  method  chosen  to  solve  each 
subset  of  equations,  that  of  Gaussian  elimination  (see  section  3.3.9  for 
rationale  and  form  of  implementation).  For  present  purposes,  the  impor¬ 
tant  feature  of  Gaussian  elimination  is  that  the  equations  are  formulated 
in  such  a  way  as  to  eliminate  the  derivative  in  eq  (12).  Then  eq  (14) 
takes  the  form 


r  =  "N  3 

s  s  s 


and  (13)  is  the  only  normal ization  condition.  If  Ts  is  large  ana  Bs 
is  small  compared  to  some  weighting  of  the  a's  in  (12),  then  the  total 
population  of  the  set  may  increase  by  orders  of  magnitude.  If  Bs  is 
large  and  Ts  small,  the  total  population  may  decrease  by  orders  of  mag¬ 
nitude.  If,  as  is  always  the  case  when  populations  approach  steady  state, 
both  Ts  and  Bs  are  small,  then  the  equation  set  approaches  indeter¬ 
minacy  and  Ns  approacnes  zero  for  al  1  s. 

To  avoid  this  problem,  Ts  is  compared  to  the  bound-Douna  input 
terms  for  that  state  of  the  set  with  the  largest  population,  s', 


Tr ,  <  z  E 


°ks'Nk 


and  3S  is  compared  to  bound-bound  output  terms. 
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If  either  test  fails  then  the  state  s'  is  removed  from  the  set,  held 
constant,  and  becomes  a  source  and  sink  of  electrons  for  the  remaining 
members  of  the  set,  so  that  they  will  have  appropriate  values  of  Ts, 

Bs. 

The  current  value  of  e  is  10“ 6  but  the  precise  value  is  not 

crucial . 

3.3.7  Sorting  Algorithm. 

The  sorting  algorithm  currently  in  use  is  applied  separately  to 
DE  and  QS,  potentially  breaking  each  of  these  major  groupings  into  smaller 
sets  of  equations. 

The  first  requirement  is  that  all  members  of  a  set  be  contiguous 
in  index,  which  is  equivalent  to  being  contiguous  in  energy  since  the 
indexing  system  for  eigenstates  runs  from  low  energy  states  to  high.  This 
condition  produces  some  unnecessary  sets  occasionally  but  prevents  the 
possibility  of,  say,  including  the  ground  states  of  two  ions  in  a  single 
set  relegating  the  excited  states  connecting  then  to  another  set.  Should 
such  a  situation  occur,  the  intermediate  states  may  not  be  able  to  come  to 
the  correct  configuration  since  they  will  be  unable  to  alter  the  ratio 
between  the  two  ground  states. 

The  second  condition  for  membership  of  state  i  in  a  set  is 
applied  only  to  QS  equations.  It  requires  that  for  some  member  j  of  the 
set, 

( aij  +  aj  i )  At  >  G  (21 ) 

This  test  is  similar  to  that  of  (13)  and  the  same  value  of  G  is 
currently  used.  The  difference  is  that  it  insists  that  state  i  meet  the  QS 
test  for  some  individual  state  j  in  the  set,  rather  than  that  i  meet  the 


test  globally  through  the  sum  of  ct-jj  over  all  j.  It  prevents  the  possi¬ 
bility  of  i  being  QS  with  respect  to  some  member  outside  the  set  but  not 
to  members  of  the  set.  An  example  is  the  upper  state  of  charge  state  Z 
being  QS  because  it  is  coupled  tightly  to  the  ground  state  of  Z+l,  which 
may  be  handled  as  a  differential  equation,  and  upper  state  of  Z  not  cou¬ 
pled  as  tightly  to  any  other  state. 

The  final  criterian  applies  to  both  DE  and  QS.  It  is  that  for 
state  i  to  be  included  in  the  set,  it  must  be  true  that  for  some  state  j 
in  the  set  either 

ay  aik  for  all  k  (where  k  can  be  DE  or  QS) 
or 

i,  i  _>  for  all  k. 

That  is,  the  transition  rate  from  i  to  some  member  of  the  set 
must  oe  the  largest  transition  rate  out  of  i,  or  the  transition  rate  from 
some  member  of  the  set  to  i  must  oe  its  largest  transition  rate. 

These  criteria  sound  quite  restrictive  and  may,  in  fact,  oe 
overaone.  However,  because  of  the  tendency  of  eigenstates  to  transfer 

most  rapidly  to  an  adjacent  state,  the  results  of  application  are  not  as 
restrictive  in  practice  as  might  seem  to  be  the  case.  Breakpoints  in  sets 
typically  occur  where  one  wants  them,  at  states  with  dipole  forbidden 
transitions  to  neighbors  and  between  states  for  which  data  exists  and  the 
artificial  uppermost  state  of  an  ion. 

When  electron  density  is  very  low  a  situation  occurs  where  the 
largest  rate  for  many  excited  states  is  spontaneous  emission  to  the  ground 
state.  Then  most  sets  consist  of  a  single  state  which  is  advanced  by  the 
time  splitting  method  described  in  3. 3. 6. 2.  No  proDlem  ensues,  not  even  a 
noticeable  penalty  in  CP  time,  for  in  this  situation  the  equations  have 
benign  behavior  and  the  time  splitting  technique  is  faster  than  Gauss 
elimination,  the  method  used  to  redistribute  population  within  a  set. 


3.3.8 


Freezing  Problem. 


The  elements  of  the  method  for  solution  can  be  summarized  as 
follows.  In  order  to  avoid  stiff  differential  equation  problems  the  equa¬ 
tions  for  eigenstate  population  were  broken  into  two  sets,  one  set  of 
stiff  equations  (QS)  for  which  solution  methods  exist  and  another  set  of 
limp  equations  (OE)  for  which  other  methods  of  solution  exist.  However, 
the  best  methods  we  have  found  are  imperfect.  To  counter  these  imperfec¬ 
tions  the  equations  may  be  further  sorted  into  smaller  sets. 

Now  what  we  would  like  to  do  in  the  interest  of  simplicity  and 
cleanliness  of  logic  is  redistribute  eigenstate  populations  within  each  of 
the  final  sets  as  influenced  by  initial  conditions  of  all  other  sets;  then 
renormalize  the  population  of  each  set  to  its  initial  value  to  prevent  an 
instability  due  to  overcorrection;  and  finally  to  adjust  the  relative  set 
populations  with  the  time  splitting  technique.  This  approach  works  well 
for  almost  all  situations  but  fails  in  situations  similar  to  the  follow- 
i  ng . 


Suppose  ion  Z+l  is  recombining  to  ion  Z.  Also  suppose  equation 
set  -i  consists  of  the  lowest  states  of  Z+l  and  perhaps  the  highest  states 

of  Z  and  is  highly  populated.  Set  A  is  rapidly  feeding  electrons  into  set 
B,  the  intermediate  states  of  Z.  Set  B  has  a  small  population  despite  the 
electrons  received  from  A  because  it  rapidly  feeds  them  downward  to  set  C, 
the  lower  lying  states  of  Z.  Set  C  is  also  highly  populated.  Further, 
suppose  set  B  is  currently  some  how  underpopulated.  Situations  similar  to 
this  are  common  for  either  recombination  or  ionization.  They  present  a 
choice  between:  (i)  a  stiff  equation  problem  where  the  time  step  is  con¬ 
trolled  by  insisting  that  the  population  of  the  low  population  set  B  not 
change  too  much  in  a  cycle  or  (ii)  a  stability  problem  where  the  popula¬ 
tion  of  set  B  fluctuates  wildly  wnile  sloshing  electrons  erratically 
between  sets  A  and  C. 


Application  of  our  desired  method  can,  under  exceptionally  un¬ 
lucky  circumstances,  proceed  as  follows  if  the  time  step  is  chosen  to  be 
appropriate  to  the  A  to  C  transfer,  as  one  would  wish,  rather  than  the 
much  smaller  time  step  appropriate  for  B  from  A  or  B  to  C. 

First,  the  sets  are  solved  independently  using  initial  condi¬ 
tions.  This  results  in  an  appropriately  modest  depletion  of  set  A,  an 
enormous  relative  net  increase  in  population  of  B  (since  it  was  assumed 
underpopulated  -  it  gains  more  from  A  than  it  loses  to  C  and  both  the  gain 
and  loss  are  very  large  compared  to  the  population  of  B  because  of  the 
large  time  step),  and  a  negligible  increase  in  the  population  of  C  (essen¬ 
tially  equal  to  the  too  low  initial  population  of  B). 

The  next  step  of  the  method  renormalizes  all  three  sets  back  to 
their  original  population,  nipping  this  incipient  instability  in  the  bud. 
The  result  is  red i str ibuted  eigenstate  populations  within  each  set  but 
initial  total  set  populations. 

Now  the  method  looks  at  the  effective  set-to-set  transition 
rates  and  picks  out  the  maximum,  say  it  is  B  to  C.  New  B  and  C  popula¬ 
tions  are  then  found  from  (17);  if  the  time  step  is  as  large  as  assumea  it 
is  likely  that  B  and  C  will  come  to  the  steady  state  ratio  they  would  have 
if  state  A  weren't  present,  which  keeps  set  B  population  below  its  correct 
value,  and,  since  B  started  with  too  low  a  value,  not  enough  electrons 
have  oeen  transferred  to  C.  Sets  B  and  C  are  then  combined  into  a  set 
(B,C)  with  relative  8  to  C  population  frozen.  New  rates  are  calculated 
between  A  and  the  combined  sets  (B,C)  and  finally  the  two  sets  A  and  (B,C) 
solved  with  (17) . 

To  the  degree  that  electrons  are  properly  distributed  in  all 
three  sets  the  result  is  that  A  is  correctly  depleted,  C  is  correctly 
enhanced  except  for  the  incorrect  sharing  of  electrons  between  B  and  C 


(minor  because  B  has  low  population),  and  B  is  left  with  a  population 
frozen  at  too  low  a  value  with  no  way  to  build  it  up  to  the  proper  value. 
In  fact,  because  the  population  of  8  is  kept  too  small,  the  electrons  in  B 
will  be  incorrectly  distributed  and,  to  some  extent  this  is  also  true  for 
A  and  C  so  neither  the  set  populations,  the  electron  distributions  nor  the 
rates  will  be  quite  right. 

There  are  two  mathematically  sound  ways  to  cure  this  sort  of 
problem,  (i)  reduce  the  time  step  to  a  value  which  produces  a  small  rela¬ 
tive  change  in  set  B  population  from  both  A  to  B  and  B  to  C  transitions 
and  (ii)  elaborate  (17)  to  handle  three  sets  simultaneously.  The  former 
is  unacceptable  because  of  CP  time  restrictions  and  the  latter  has  not 
been  tried  because  there  is  no  guarantee  that  a  more  elaborate  4  or  5  set 
problem  would  not  occur. 

What  has  been  done  instead  is  an  empirical  reduction  in  the 
degree  of  renormal ization  of  small  population  sets  after  the  initial  solu¬ 
tion.  The  problem  only  occurs  for  sets  with  small  population  and  would 
not  exist  if  sucn  sets  were  not  renormalized  after  the  initial  solution 
for  electron  distribution  within  the  set.  Therefore  the  renormalization 
procedure  is  modified  to  the  following: 

(1)  The  renormal ization  factor,  R,  is  calculated  for  the  set 

(2)  If  the  set  population  is  larger  than  some  fraction,  f,  of 
the  total  element  (all  Zs)  population,  renormalize.  The 
present  value  of  f  is  1/3. 

(3)  If  the  set  population  is  less  than  f  then  replace  R  by  the 
function 

R'  =  Ra+l-a  (22) 


and  renormalize  using  R1 2 3. 


The  function  a  is  given  by 


a  =  exp(|2-R-l/R|  )  (23) 

From  (23)  one  sees  that  for  R  »  1  or  R  «  1  the  function  a  ap¬ 
proaches  zero  and  R1  therefore  approaches  1.  In  this  case  almost  no 
renormal i zation  takes  place,  the  modified  value  of  set  population  is 
accepted,  allowing  low  population  sets  to  approach  their  proper  popula¬ 
tion.  On  the  other  hand,  when  R‘  is  near  unity,  a  approaches  unity  faster 
than  R  and  renormal ization  is  accepted,  preventing  oscillation  about  the 
correct  value. 

The  approach  is  inelegant  but  yields  good  answers  without  signi¬ 
ficant  penalty  in  wasted  CP  time. 


3.3.9  Gauss  Elimination. 

3.3.9. 1  Quasi-Steady  Equations.  In  the  quasi-steady  approximation  the 
derivative  in  (12)  is  small  compared  to  other  terms  and  is  neglected.  The 
set  of  equations  reduces  to  a  set  of  linear  algebraic  equations  for  1  _<  s 
■:  m, 


N  (S  +B  )  =  T  +  I 
sws,m  s'  s  s,m 


(24) 


where  for  convenience  we  have  defined  the  bound-Douna  output  term 
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as,s  =  0 


(25) 
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and  the  bound-bound  input  term 
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(26) 


There  are  at  least  three  standard  approaches  to  solution  of  the 
set  (24):  determinants,  iteration,  and  Gauss  elimination.  Both  determi¬ 
nants  and  Gauss  elimination  are  quite  sensitive  to  round  off  errors  so  our 
initial  choice  was  iteration  despite  the  fact  that  it  is  reputed  to  be  the 
slowest  of  the  three  methods.  The  allegation  was  amply  verified  by  our 
experience  and  eventually  we  were  driven  to  Gauss  elimination. 

The  advantage  of  Gauss  elimination  is  that  it  is  faster  than 
iteration.  Its  advantage  relative  to  the  method  of  determinants  is  that 
it  was  easier  to  analyze  its  sensitivity  to  rouna  off  problems.  In  fact, 
we  were  able  to  devise  a  modification  of  the  technique  which  completely 
eliminates  this  problem  for  equation  systems  like  (24). 
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To  see  how  to  do  this  consider  the  mechanics  of  Gauss  elimina¬ 
tion.  Initially  one  has  the  m  equations  (24)  in  m  unknowns,  Ns .  The 
set  is  reduced  to  m-1  equations  in  m-1  unknowns  by  substituting  the  mth 
equation  into  all  the  others,  then  to  m-2  equations  in  m-2  unknowns  by 
substituting  the  m-i  equation  into  the  remainder,  etc.,  until  only  one 
equation  in  one  unknown  (N j)  remains.  This  can  be  solved.  Then  is 
substituted  to  find  N2  etc.,  until  all  Ns  are  found. 

Given  that  all  terms  of  (24)  are  nonnegative  as  is  our  case,  the 
only  way  to  produce  a  serious  round  off  error  is  by  some  step  which 
requires  a  subtraction.  This  occurs  when  substituting  the  mth  equation 
into  any  equation  s  as  follows. 
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32 
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However,  further  study  of  (33)  shows  that  the  subtraction  can  be 
eliminated.  Look  more  closely  at  the  second  and  fourth  terms  on  the  RHS 
of  (33) 

a._  -  cl,.  acm/(Sm  +B  ) 
sm  ms  sm  m,m  m 

=  a  (S  +B  -a.  )/(S  +B  ) 
sm'  m,m  m  ttis  m,m  m 

=  a._(S_  i+cl+S_  ..i  +B_-cl,.  )  /  ( S  +B  )  .  (34) 

sm'  m,s-l  ms  m,s+l  m  ins  m,m  nr 


The  only  subtraction  cancels! 


In  fact,  we  can  use  the  same  modification  of  a  as  in  (32)  to 


find 


m-1 


Ss,m-1  =  as,k 


(35) 


and  define 


Bl 


B  +  a  B  /  ( S  +B  ) 

s  sm  m  '  m,m  nr 


(36) 


and  Eq.  (29)  has  the  proper  form  with  no  subtractions.  Accordingly  no 
round  off  problem  exists  until  the  number  of  ecuations  approaches  the  big¬ 
gest  integer  which  can  be  represented  by  the  computer  word  length  ( 1 0 7  for 
the  smallest  modern  machines),  or  perhaps  the  square  of  the  biggest  inte¬ 
ger  (1014). 


Indeed,  back  substitution  of  answers  into  the  original  equa¬ 
tions,  when  m  is  around  40,  show  discrepancies  of  at  most  the  last  three 
bits  for  a  32  bit  word. 


3. 3. 9. 2  Differential  Equations.  Eq.  (12)  can  be  converted  to  a  simple 
first  order  difference  scheme  by  substituting 

dN$/dt  =  (N  (t+At)  -  N s ( t ) ) / At  (37) 

Then  (12)  can  be  written 


Vt+lt>  <s;,m  *  8s> 


t:  *  [■ 
s  s  ,m 


which  is  identical  in  form  to  (24),  if  we  define 


ak,s  5  At  “ks 


S'  =  V 
s,m  s>l< 


3'  i  AtS 
s  s 


r  =  N$(t)  +  AtTs 


V  Vks 


Having  done  this,  the  same  Gauss  elimination  routine  as  was  used 
for  quasi-steady  equations  can  be  used  to  obtain  a  first  order  implicit 
solution  for  differential  equations. 


"  •  *  .  *  * .  .  •  .  •  «.**.  • 


SECTION  4 

RADIATION  SPECTRUM  AND  LINE  SHAPE 


4.1  NOMENCLATURE. 

One  speaks  of  of  "lines"  meaning  spectral  emission  or  absorption 
features  as  well  as  energy  quanta  absorbed  from  the  spectrum  or  emitted  to 
the  spectrum  by  an  ion.  When  translated  into  Fortran  such  usage  can  be 
very  confusing.  Therefore  in  the  following,  the  radiation  spectrun  repre¬ 
sented  by  a  set  of  numbers  denoting  frequencies,  energy  flux,  etc.  has  no 
lines.  It  has  "features."  A  "line"  will  be  associated  with  an  atom  or 
ion,  it  refers  to  energy  quanta  absorbed  from  the  radiation  field  or  emit¬ 
ted  to  that  field  associated  with  transitions  among  electronic  eigen¬ 
states. 


4.2  SPECTRAL  FEATURES. 

The  convention  adopted  is  to  represent  the  radiation  field  as 
being  composed  of  a  number  of  spectral  "features,"  hopefully  less  than  100 
of  them  in  a  real  problem,  although  1000  are  allowed.  Each  feature  is 
defined  by  four  parameters;  its  central  energy,  v^(eV),  its  central,  or 
core  intensity,  Uc(erg/cm2/sec/eV)  or  (erg/cm2/sec/eV/sr) ,  its  core  half 
width,  Wf(eV),  and  a  parameter  denoting  the  strength  of  its  tail, 
Ul(erg-eV/cm2/sec)  or  (erg-eV/cm2/sec/sr) .  The  feature  intensity  is 
calculated  from  parameters  as 


(  UT/(v-vf)2 


if  |  v-vw|  >  Wf 

In  addition  to  these  four  parameters  there  are  three  rules  of 
logic:  (1)  each  feature  is  associated  with  a  "level",  level  1  feature 
parameters  are  stored  first,  followed  by  level  2,  etc.  Levels  are  distin¬ 
guished  on  the  basis  of  core  width;  roughly  speaking,  core  width  should 
become  narrower  by  at  least  about  a  factor  of  four  when  trie  level  in¬ 
creases.  That  is,  higher  levels  represent  narrow  features  superimposed  on 
broader  ones,  (2)  within  a  given  level,  features  are  stored  sequentially 
on  the  basis  of  central  energy,  \y,  (3)  level  1  represents  the  continuum; 
its  features  have  Uy  =  0.  There  are  9  levels  allowed.  There  is  no 
limit  on  the  numoer  of  features  alloweu  in  each  level,  only  on  the  total 
number  of  features.  The  logic  at  present  requires  at  least  one  feature  to 
oe  in  level  i,  althougn  its  intensity  may  be  zero. 

The  system  is  believed  capable  of  representing  a  radiation  field 
to  any  desired  degree  of  accuracy  needed  for  energy  transport  and  the 
associated  problem  of  eigenstate  population  as  .veil  as  output  spectrum  for 
systems  analysis  and  data  comparison. 


LINES. 


Atomic  (ionic)  lines  are  also  modeled  with  four  parameters, 
'■epresenting  a  Doppler  oroadened  line  core  with  Lorentzian  wings.  The 
parameters  are  central  energy,  vj_(eV),  Lorentzian  central  cross  section 
q_(cm2),  Lorentzian  half  width,  W(_(eV),  and  the  Doppler  width  to 
L.orentzian  width  ratio,  x. 


W  =  U  /7 


a  =  a.  /  /l  +X  2/  it 

c  L 


2.490606  +  0.539095  x1*  65668 


1.74623  x1*  08419 


,  x  _<  4 
,  x  >  4 


j(  v)  = 


ac  exp  [-(  v- 2/W 2  ] 


l+(  v-  )  2/W^  2 


.  ini 


- <  X 

u  —  C 

WL 


lriL 


This  representation  gives  the  correct  central  cross  section  ana 
the  correct  cross  section  in  the  far  wing  always.  For  x  less  than  about  3 
it  has  a  maximum  error,  as  compared  to  the  more  correct  Voight  profile,  of 
under  30  percent,  occurring  at  the  transition  point,  xc.  When  x  >  3  the 
error,  an  underest imate  of  cross  section,  becomes  worse,  but  when  x  >  3, 
the  wing  portion  is  less  than  10” 4  of  the  center  so  this  error  cannot 
seriously  affect  any  overlap  integral,  thus  cannot  seriously  affect  any 
rate. 

4.4  LINE  WIDTHS. 

The  Doppler  width  of  a  line  is  given  by 


WQ  =  v^/c  /2  0../M 


,  V,  /.V.V.  <•. 

.•j'V ->  v. 


where  c  is  light  speed  in  vacuo,  9i  is  heavy  ion  temperature,  and  M  is 
heavy  ion  mass. 

The  lorentzian  width  is  the  sum  of  widths  of  the  two  eigenstates 
involved  in  the  transition.  For  each  of  these  states,  width  is 

WL  =  Wr  (44) 

where  r  is  the  effective  lifetime  of  the  state  of  ion  Z  against  all  forms 
of  decay,  plus  thermal  electron  induced  Stark  broadening,  rs. 

r  ■  rs  *  %  *  rEB  *  rEF  *  rER  *  rPF  *  rPR  *  rPB  <45> 

where  rsp  is  spontaneous  emission  rate  to  all  lower  states  of  Z,  r^B  1S 
electron  collisional  oound-bound  rate  to  all  states  of  Z,  rpp  is  elec¬ 
tron  col  1 isional  ionization  rate  to  all  states  of  Z+l,  rpR  is  electron 
three-body  recombination  rate  to  all  states  of  Z-l,  Tpp  is  pnoto  ioniza¬ 
tion  rate  to  all  states  of  Z+l,  Tpr  is  radiative  recombination  rate  to 
all  states  of  Z-l,  and  rpg  is  bound-Dound  photo  induced  transition  rate 
to  all  states  of  Z. 

These  rates  are  calculated  in  the  order  shown  in  (45),  from  left 
to  right,  because  the  rates  become  increasingly  sensitive  to  state  width 
from  left  to  right.  In  principle  one  should  iterate  (45)  to  obtain  a  self 
consistent  state  width.  In  practice,  the  bound-bound  photo  transition 
rate  is  the  most  sensitive  to  width;  if  rps  changes  r  by  as  much  as  25%, 
that  one  contribution  is  iterated. 

The  method  for  calculating  each  of  these  rates  is  described  in 


SECTION  5 
RATES 


In  the  following  discussion  the  convention  is  adopted  that  the 
parenthetical  expression  (a,b)  means  "from  a  to  b." 

5.1  STARK  WIDTH. 

The  net  effect  of  Stark  shifts  due  to  bombarding  thermal  elec¬ 
trons  is  to  broaden  eigenstate  energy.  This  process  does  not  cause  tran¬ 
sitions  between  states  but,  because  its  effect  on  line  width  is  similar, 
it  is  discussed  here.  To  a  first  approximation  Stark  broadening  is  not 
dependent  on  state  width,  so  it  is  appropriate  to  calculate  it  oefore 
other  contr ibutions  to  state  width  have  been  evaluated. 

A  formula  for  thermal  electron  Stark  broadening,  rs,  is  given 
in  eq  (526)  of  Ref.  2.  This  formula  can  be  fit  within  a  few  percent  over 
the  interesting  range  of  electron  temperature,  e,  by 

rs  =  W$Ne//0  [1  +  r/(As+r)  +  C$ 0Ps  ]  (46) 

where  r  is  the  ratio  of  e  to  the  Rydberg  (13.6  eV)  and  W$,  A  ,  Cs,  and  p$ 
are  parameters  adjusted  to  give  a  best  fit  to  Griem's  formula  for  each 
eigenstate.  These  parameters  are  derived  at  the  time  eigenstate  data  is 
collected  and  treated  as  eigenstate  data  (see  section  6). 


5.2 


SPONTANEOUS  EMISSION. 


Spontaneous  emission  is  also  insensitive  to  state  width  and  is 
the  second  contribution  to  state  width,  the  first  actual  transition  rate, 
calculated.  From  Ref.  3.,  p.  58,  the  spontaneous  radiative  transition 
rate  from  upper  state  k  to  lower  state  i  is 

rsp(k,*)  =  (8iT2e2v2)/(mc3)(gji/gk)f  &  (47) 

in  c.g.s.  units.  In  (47),  rsp  is  in  units  of  sec'1,  e  (e.s.u.)  is 
electron  charge,  m(g)  is  electron  mass,  c  (cm/sec)  is  lightspeed,  the  g's 
are  state  degeneracies  and  f^  is  the  oscillator  strength  from  state  i 
to  state  k.  In  more  convenient  units  (47)  becomes 

■VM)  ■  V£k-£t)!  <9A)f«k  1481 

wnere  the  E's  (eV)  are  state  energies  and 

Af  =  4.33927xI07  eV-2  sec'  1  (49) 

5.3  ELECTRON  BOUND-BOUND  COLLISIONS. 

Burgess  (Ref.  4)  gives  a  formula  for  electron  collision  cross 
section  from  lower  state  i  to  upper  state  k  as 

^(C«2)  ■  ^?ftk/lE(Ek-Et)!  ,  £  >  £*"£;  (50) 

where  A  =  1.28x10"  15  cm2,  R  is  the  Rydberg,  g  is  the  Gaunt  factor,  f  i 

y  kx 

the  oscillator  strength,  E  is  electron  energy  and  E^,  E  are  the  state 
energies 


To  convert  (50)  into  a  rate  for  electrons  at  temperature  e  one 
multiplies  by  electron  velocity  then  averages  over  a  Boltzmann  distribu¬ 
tion  to  obtain 

rEB(i,k)  =  CEC  f)lk/[(01/2(Ek-Ejl)]  exp(-(Ek-EJt)/0)  (51) 
CEC  =  6x10“ 6  eV3/2  sec"1 


where  g  was  taken  to  be  0.4  and  9  and  E  are  in  units  of  eV. 

By  detailed  balance,  the  downward  rate  from  k  to  t  is 
rEB(k,<0  =  CECf  ilk9 ^9k 01/2  (Ek"Et^ 

5.4  ELECTRON  BOUND-FREE  COLLISIONS  AND  THREE  BODY  RECOMBINATION. 

In  Ref.  6.,  Lotz  gives  a  general  formula  along  with  necessary 
parameters  which  gives  a  good  fit  to  measured  electron  collisional  ioniza¬ 
tion  cross  sections  and  rate  coefficients  for  any  atom  or  ion  in  the  lower 
portion  of  the  periodic  table.  His  fit  is  based  on  the  ratio 

x  =  aE/9  (53) 

where  &E  is  change  in  energy,  ground  state  Z  to  ground  state  Z+l.  He  sug¬ 
gests  that  ionization  rates  from  excited  states  be  approximated  by  simply 
substituting  the  appropriate  value  of  aE.  We  have  adopted  this  suggestion 
and  generalized  it  to  include  excited  states  of  Z+l  as  well. 


For  our  application,  his  formula  can  be  well  fit  to  the  simpler 


formul  a 


rEF(k,p)  =  (AL/e3/2)NeE  j(x)/x 


L 


(54) 


where 


X  =  (Ep+Va-Ek)/ 0  ,  (55) 

9  is  electron  temperature,  Ne  electron  density,  Ej  the  first  order  expo¬ 
nential  integral,  Ep  and  E|<  are  state  energies  in  Z+l  and  Z  respec¬ 
tively.  Va  is  the  effective  ionization  energy  of  Z,  and  AE  and  PE 
are  parameters  fit  to  Lotz's  formula  for  ion  Z.  These  parameters  are 
derived  by  comparing  the  sun  of  (54)  over  p  to  the  Lotz  formula  for  state 
k  of  Z. 


It  is  easy  to  criticize  this  procedure  out  criticism  can  oe 
readily  countered  by  comparison  of  the  formula  to  data  except  for  the  gen 
eralization  to  excited  states  of  Z+l  (for  which  there  is  no  data).  That 
criticism  must  be  countered  by  a  challenge  to  find  a  better  formula, 
coupled  with  the  observation  that  ionization  rates  to  excited  states  of 
Z+l  are  almost  always  too  low  to  control  state  populations. 

The  reverse  reaction,  three-body  recombination,  can  be  derived 
from  rEE  by  use  of  the  principle  of  detailed  balance,  which  yields 

rEB(p,k)  =  rNegk  exp(x)  Vrqe93/2gpl  r£F(k,p)  (56) 

Q  =  6.037xl021  eV"3/2  cm3 

where  g^  and  gp  are  statistical  weights  and  Qg  is  the  partition  function 
of  a  free  electron  at  i-eV  temperature. 


5.4.1 


Contribution  from  Eigenstates  Above  the  Continuun. 


The  quantity  Va  in  (55)  is  the  ionization  potential  depressed 
by  polarization  effects  of  the  plasma.  In  all  our  formulae,  states  which 
lie  above  Va  are  neglected.  But  these  eigenstates  still  exist  and  the 
oscillator  strength  to  them  still  exists  so  they  may  possibly  have  some 
influence.  To  evaluate  the  magnitude  of  influence  of  such  states  we  con¬ 
structed  and  implemented  a  simple  model.  It  goes  as  follows.  The  infinite 
number  of  eigenstates  which  lie  between  Va  and  the  vacuum  ionization 
potential,  V,  all  lie  above  the  ground  state  energy  of  ion  state  Z+l,  in 
most  cases  extremely  close  to  this  ground  state.  So  assume  any  bound- 
bound  electron  collision  to  one  of  these  states  to  be  instantly  followed 
by  another  collision  which  frees  the  electron,  resulting  in  a  ground  state 
1+1  ion. 


To  evaluate  this  model  we  need  to  sum  (51)  over  all  states  above 

3iTI  and  add  the  resulting  rate  to  r^p(k,l)  (from  Eq.  54),  that  is,  to 

the  electron  col  1 isional  ionization  rate  to  the  ground  state  of  Z+l. 

If  all  states  are  assumed  hydrogenic,  Kramers  rule  (65)  is 

adopted,  and  in  the  exponential,  upper  state  energy  is  taken  to  be  7a  to 
produce  a  simple  overestimate  of  the  effect,  one  finds 

rt,  l  =  2  CccVtVe1/2)  exPr-(Va-Ez)/eU3I  (57) 


wnere 


I  =  V  k  5/  ( k  2-  jz2)  u 
5+1 

71 


44 


A  reasonably  good  approximation  to  I  may  be  found  by  defining 


a  *  Sm  • 


,  b  i  (S  +1 )  2-  l2  ,  c  =  (S  +1.7)2-  £s 


I  =  S|/a4  +  (Sn+l)5/b4  +  1/2  [U4/3)/c3  +  i2/c2  +  1/c] 

The  term  (57)  is  added  to  the  electron  bound-free  rate  when  p  = 
1.  Its  most  important  effect  should  occur  at  high  electron  density  when 
Sm  becomes  low. 

To  date  no  case  where  the  additional  rate  (57)  is  important  has 
been  found.  This  has  discouraged  us  from  seeking  corresponding  effects 
for  other  processes.  Further,  unless  we  are  able  to  find  some  case  of 
potential  interest  for  which  the  effect  matters,  we  will  delete  it  from 
the  code  to  save  time,  storage,  and  eliminate  unnecessary  complication. 

5.5  RADIATIVE  IONIZATION  AND  RECOMBINATION. 

5.5.1  Radiative  Recombination. 

The  radiative  recombination  coefficient  from  state  k  of  ion  Z  to 
state  l  of  ion  Z-l  can  be  obtained  from  the  inverse  photoionization  cross 
section  through  the  principle  of  detailed  balance,  but  not  as  simply  as 
for  collision  processes. 

Let  aU,Z-l;k,Z; v)  e  a  be  the  photoionization  cross  section  at 
frequency  •>  from  i  to  k.  a(k,Z;  £,Z-1 )  5  a  be  the  radiative  recombination 
coefficient  from  k  to  i. 


Then  in  general,  for  these  processes  only. 


oo 

N  =  cfleN(Z,k)  -  N  /  oU(  v)d  v  (59) 

Vo 

where  Ng  is  electron  density,  N(Z,k)  is  the  density  of  Z  ions  in  state  k, 
N  is  density  of  Z-l  ions  in  state  i,  is  the  minimum  frequency  which 
can  ionize  from  n  to  k  and  U( v)  is  the  photon  omnidirectional  spectral 
flux  (photons/cm2/sec/hz) . 

To  apply  detailed  balance,  assume  (59)  to  be  written  for  equi- 
librium  conditions.  Then  N^must  be  zero  and  ll(  v)  must  be  Planckian, 

U(  \>)  =  P(e)  =  Cp  e2/(exp{  t) -1)  (60) 


where 


t  =  e/9 

Cp  =  (8*eJ)/(h3c2) 


e(eV)  is  photon  energy,  e(eV)  is  temperature,  e0  converts  from  eV  to  ergs, 
n  is  Planck's  constant  and  c  is  the  speea  of  lignt.  Substitute  (60)  for 
U(v)  in  (59)  and  set  N*  =  0  to  obtain 


NeN(Z,k) 


r  -f{  e)  d  e 

^0 


(6i ; 


The  population  ratio  in  front  of  the  integral  can  be  found  for 
aqu  i  1  i  or  i  urn  conditions  tnernooynain  ical  i  y ,  Dy  taxing  the  ratio  of  partition 

functions 


i 


8 


NeN(Z,k)  Q(e)  Qk 


=  h^exp^/fg^UumkT)3/2]  (62) 


In  our  units  this  is 


yirfer  Se«P(*)/t9k<Je«3/2l 


where  Qg  =  6. 03 704x10 21  cm- 3  eV-3/2  and 


(E,+I,_,-E.)/e 


*  =  ukT1z-rc£ 


is  the  ratio  of  energy  required  to  ionize  state  i  then  bring  the  ion  to 
state  k,  to  the  temperature. 


Now  we  have  the  multiplier  of  the  integral  in  (61)  in  terms  of 
known  quantities  and  must  evaluate  the  integral.  To  do  that  requires  an 
expression  for  the  photoionization  cross  section. 


3egin  with  Kramer's  formula  for  ionization  of  hydrogen  from 
principle  quantun  number  n 


<7k(v)  «  (64AiZ4e10)/(3/3  ch6)g/(n  V) 


wnere  n  is  electron  mass,  e  is  electron  charge  and  g  is  the  Gaunt  factor. 


Somehow  (65)  must  be  generalized  to  apply  to  ions  other  than 
nytfro gen  and  to  be  specific  with  respect  to  the  final  state  of  the  ion  Z. 


B— PI 
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We  chose  to  generalize  (65)  beyond  hydrogen  by  eliminating  as 
much  specific  Z  and  n  dependence  as  possible,  replacing  it  with  energy 
dependence,  through  the  hydrogenic  formula  for  ionization  potential 

Iz  =  (2ir2mZ2e4)/(h2n2)  (66) 

Substitution  of  (66)  into  (65)  to  eliminate  Z  dependence  and 
conversion  to  our  units  results  in 

cr(  e)  -  CHI2/(ne3)  (67) 


where 


CH  =  (16e  2hg ) / ( 3  73  mce0)  =  3.2275“  17  cm2  eV 

if  the  Gaunt  factor,  "g,  is  taken  equal  to  0.3. 

To  use  (67)  for  specific  states,  the  simplest  method  is  to  just 
replace  I  by 

*z  =  Ek  *  1  -  E*  (68) 

where  I  is  ionization  potential  of  the  ground  state  of  Z-l. 

Substitution  of  (68)  into  (67)  followed  by  substitution  of  (60), 
(63),  and  (67)  into  (61)  leads  to 

a  =  (p(  £,k)CHCp91/2g!i)/(Qengk)F(x) 


F(x)  =  x^e  /  dt/[t(e  -1)  j  ,  xQ  =  I^/e 


and  the  constant  p(  z,k)  will  be  unity  unless  data  on  a  specific  transition 
indicates  otherwise. 


An  approximation  for  F(x)  good  to  better  than  1%  everywhere  is 


x2ex  [1/x  +  ( znx) /2  -  x(l+x/2)/12  -  .62772  ]  ,  x  <  0.25 
x2exE  j(x)  [x2+a  ^+82]/ [x^b  1x+b2]  ,  0.25  _<  x  <  4.3 


F(x)  = 


x  2exE  x(  x)  , 


4.3  £  x  <  60 
60  <  x 


with  aj  =  2.9823,  a2  =  2.4416,  b  t  =  3.5329,  d2  =  0.38094  ana  E  x( x)  is  the 
first  oraer  exponential  integral. 

5.5.2  Photoionization. 

The  rate  of  photoionization  from  state  z  of  ion  Z  to  state  k  of 
ion  Z+l  is  given  by 


r  z  k  =  '  y 

v0 


where  a  is  given  by  (67)  and  U( v)  is  the  omnidirectional  spectral  flux 
described  in  Section  4.2,  eq.  (40). 


!to 

$ 


To  bring  (40)  into  proper  units,  (71)  becomes 

rt,k  =  CHe‘o1  V  vfl  IZ/n  /  U(£)/e3de 
f  f 


where  the  subscript  f  denotes  a  spectral  feature. 

The  representation  of  U  given  in  (40)  requires  three  types  of 
integrals  to  be  carried  out  in  (72)  corresponding  to  a  rising  tail  of  U  if 
e  <  a  constant  U  if  -  Wp  <  e  <  \>p+Wp  and  a  falling  tail  if 

e  >  v^+Wp.  If  the  threshold  energy  for  ionization,  e-,  lies  below  the 
feature  core,  all  three  types  of  integral  are  present  etc. 


The  simplest  type  is  that  of  constant  U,  applicable  to  the 


core.  Then 


f  J de/e3  =  Uc/2(1/B2-1/T2) 


where  B  =  max  (  -  Wp,  Gy),  T  =  vp  +  Wp 

The  two  types  of  feature  tail  integrals  can  be  handled  with  a 
single  equation. 

If  below  the  core,  (rising  tail), 

X  =  t-r 


y  =  vp/(vf-Wp) 


If  above  the  core,  (falling  tail), 


x  max  [  vp/  Gy ,  Vy/ (  vp^*Wy) 


y  =  0 


/  Ude/  e3  =  (UT/vf4)  {3  tn  [A± Ill]  +  2  (x-y) 

+  (x2-y2)/2  +  (y-x)/[(l-y)(l-x) ]  }  (74) 

S.6  BOUND-BOUND  RADIATIVE  RATES. 

To  calculate  the  transition  rate  from  lower  eigenstate,  l,  to 
upper  eigenstate,  k,  of  ion  Z  due  to  photon  induced  transitions,  the  inte 
gral 

=  /  oUdE  (75) 

must  be  evaluated  using  the  line  model  (41)  for  a  and  the  feature  model 

(40)  for  U. 

Given  rik,  the  inverse  rate  may  be  obtained  from  detailed 
balance  and  is 

r'<  l  '  <9A>rtk  (76» 


The  integral  (75)  may  be  sensitive  to  the  line  width,  (45)  and 

the  terms,  rpg,  due  to  photo  induced  transitions  in  one  or  both  of  the 

states  i  and  k  may  dominate  line  width.  Thus  the  rate  r  ,  (and  r.  ) 

Ik  kr 

implicitly  depends  upon  itself  as  well  as  all  other  rates  to  and  from  i 
and  k.  To  partially  account  for  this  without  incurring  too  much  CPU 
penalty,  if  the  rates  r^  and  rki  increase  the  line  width  by  as  much  as 
25%,  the  rates  r^  and  r^  are  recalculated  at  the  new  line  width.  This 
will  not  cover  all  cases,  but  should  cover  all  rapid  transitions,  which 
are  the  ones  of  most  importance. 


•  •  v. 


•i  Vfr»  * 
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Evaluation  of  the  overlap  integral  (75)  is,  of  course,  the  heart 
of  the  rate  calculation.  The  line  and  feature  models  were  designed  to  be 
realistic,  yet  allow  analytic  evaluation  of  the  overlap  as  far  as  possi¬ 
ble.  The  sketch  below  illustrates  the  problem. 


■"f-pf- 


The  central  line  frequency,  may  be  located  anywhere  with 
respect  to  the  central  feature  frequency,  Vp  and  x^  may  be  larger, 
smaller,  or  identical  to  Wp. 

Eight  logically  distinct  cases  are  recognized.  Two  are  handled 
by  simple  analytic  formulae  and  the  other  six  by  various  sequences  of 
calls  to  four  subroutines.  Three  of  the  four  subroutines  evaluate  ana¬ 
lytic  formulae;  alas,  the  fourth  is  a  numerical  integration.  Perhaps 
someday  we  will  find  an  adequate  analytic  approximation  for  this  one  also. 


I«CW 


Case  1 


Line  versus  Continuum 


If  the  feature  is  in  level  1  of  the  spectral  description  it  has 
no  tail,  Uy  =  0,  and  can  be  assumed  very  broad  with  respect  to  the  line. 
Thus  in  (75)  one  can  take  U  to  be  constant  at  the  value  Uc  and  use  infi¬ 
nite  limits  for  the  integral.  In  that  case  an  exact  integral  exists  for  a 
Voight  profile  (which  is  a  better  representation  than  (45)) 

/aUdE  =  ttUc (  v|_£0)  (77) 

where  e0  converts  ergs  to  eV  and  the  symbols  are  defined  in  (41)  and  (42). 
Case  2  Line  and  Spectral  Feature  Well  Separated 


This  is  another  especially  simple  case,  and  snould  be  appropri¬ 
ate  for  most  spectral  features  for  any  given  line.  If  (Wf+W[_xc)/j  vf'v'L|  < 
0.25  then 

foUdE  =  tjW ^ a^U  f/ (  Vf-  ^l )  2  +  2a^(Uy/Wy+UcWy)/A2  (78) 


where 


A  =  1+  (v^-\^)/W^ 

The  other  cases  are:  (3)  cores  don't  overlap,  line  below  feature,  (4) 
cores  don't  overlap,  line  above  feature,  (5)  cores  partially  overlap,  line 
below  feature,  (6)  cores  partially  overlap,  line  above  feature,  (7)  fea¬ 
ture  core  contained  in  line  core,  and  (8)  line  core  contained  in  feature 
core.  Each  of  these  cases  is  handled  by  a  special  sequence  of  calls  to 
subroutines  which  evaluate  the  following  four  types  of  integrals. 


Standard  formulae  exist  for  approximating  this  integral,  one 
must  only  be  sure  to  get  the  sign  right  since  a  and  b  can  be  on  either 
side  of  and  check  for  a  and  b  on  opposite  sides  of  ^ . 

Type  4 

T4  J  e"^^  /W  /(  wf)2dv  (87) 

a 

To  the  present,  we  have  found  no  adequate  analytic  approximation 
for  (87).  It  is  integrated  numerically. 

5.7  OSCILLATOR  STRENGTH . 

5.7.1  Oscillator  Strength  for  Real  States. 

In  the  foregoing  sections,  f  ,  is  oscillator  strength  for  allow- 
ed  transitions,  and  is  input  data  for  the  calculation.  To  account  for  an 
intuition  that  nigher  multiple  moments  should  be  more  effective  for  elec¬ 
trons  than  for  photons,  we  introduce  the  arbitrary  rule  that  for  electrons 

f^  =  max  (data  for  optically  allowed  transition,  f  1 }  (38) 

where  f  is  taken  to  be  10" 3.  This  cavalier  treatment  of  optically  for¬ 
bidden  transitions  should  be  improved  someday,  but  the  above  is  how  it 
stands  at  present. 


5.7.2 


Oscillator  Strength  to  the  Artificial  State. 


In  order  to  calculate  transition  rates  between  some  real  eigen¬ 
state  with  index  i  and  the  artificial  state  with  index  a,  an  oscillator 
strength,  f  ,  must  be  calculated.  This  can  be  done  by  assuming  a  model 

ila 

for  population  distribution  among  the  component  states  of  the  artificial 

state,  summing  the  transition  rates  to  all  of  them,  and  finding  that  value 

of  f  which  yields  the  same  total  rate, 
ta  J 

The  model  assumed  is  that  eigenstate  members  of  the  artificial 
state  are  hydrogenic  and  have  a  Boltzmann  distribution  for  all  principal 
quantum  numbers  above  some  value,  b. 

The  first  step,  then,  is  to  find  b.  This  is  accomplished  by 
comparing  the  sum  of  all  spontaneous  radiative  transitions  to  lower  states 
(easy  to  approximate)  to  the  electron  coll i s i on al  excitation  rate  to  the 
next  higner  state  (also  easy  to  approximate).  The  idea  is  that  spontane¬ 
ous  radiation  is  the  only  process  tending  to  drive  state  populations  out 
of  equilibrium  and  it  drives  electrons  downward.  So  a  state  can  be  in 
equilibrium  relative  to  its  neighbors  only  if  other  processes,  in  particu¬ 
lar  those  driving  electrons  upward,  dominate  spontaneous  emission.  The 
simplest  such  process  to  consider  is  electron  collisional  excitation  from 
state  b  to  state  b+1.  It  would  be  more  elegant  and  a  little  better  to  use 
the  sum  of  all  upward  collisions  but  that  is  mathematically  cumbersome  and 
furtner,  for  hydrogenic  states  the  widest  energy  gap  to  bridge  is  to  the 
adjacent  state,  thus  that  single  transition  should  serve  as  a  good  measure 
of  the  sum  to  all  upper  states. 

Of  course,  it  is  possible  that  electron  density  or  temperature 
could  be  so  low  that  photon  induced  upward  transitions  dominate  electron 
induced  transitions  and  should  be  used  in  the  comparison.  Radiative 
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transition  rate  is  more  expensive  to  calculate  than  collisional  rate  even 
for  the  b  to  b+1  transition,  but  given  an  arbitrary  line  spectrum,  one 
cannot  assume  any  single  transition  to  be  typical,  so  the  full  sum  would 
have  to  be  approximated  somehow.  This  elaborate  process  seems  unneces¬ 
sary;  for  if  radiative  transitions  dominate,  making  our  estimate  of  b  too 
large  and  thus  our  estimate  of  f  too  small,  it  is  probable  that  radia¬ 
tion  will  be  intense  enough  to  drive  state  t  toward  equilibrium  with  state 

a  even  with  a  somewhat  low  value  of  f„,. 

ta 

If  V  is  the  vacuun  ionization  potential  of  a  hydrogenic  ion  then 
the  energy  difference  between  the  hydrogenic  state  b  and  some  lower  hydro¬ 
genic  state  l  is 

Efa  -  E  ^  =  V(b2-  l2)/b2i2  (89) 

And  Kramer's  formula  for  hydrogenic  oscillator  strength  is  aDOut 
ftb  =  2  tb  3(b  2-  t2)~  3  1 90) 


*9 


Substitution  of  (90)  and  (89)  into  (43)  and  recognition  of  the 
fact  tnat  nydrogenic  state  statistical  weights  are  proportional  to  princi¬ 
pal  quantum  number  squared  yields  the  transition  rate  due  to  spontaneous 
radiation. 


rsp(b,i)  =  2ArV2  rb  3z(b2-  i2)  1 


(91) 


To  get  the  sum  of  all  transitions  from  b,  to  states  1  through 
b-1  approximate  the  sum  of  (91)  by  an  integral.  For  b  large  compared  to 
unity  the  integral  is  approximated  by 
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s  =  y  rfb.i)  =  2A  V2/b5  fc^Benb]  (92) 

1=1  SP 

where,  if  the  normal  integration  limits  of  t  =  1/2  to  b  -  1/2  are  taken, 
a  =  .347  and  8=1.5 

A  slightly  better  fit  to  the  sun  is  obtained  for  the  most  impor¬ 
tant,  low,  range  of  b  by  taking 

a  =  0.411  ana  8  =  1.516  (93) 

Substitution  of  (89)  and  (90)  into  (51)  yields  the  electron  col- 
lisional  rate  from  b  to  b+1.  After  clearing 

rEB(b,b+l)  =  (2CEcNe/91/2V)b3(b+l)5/(2b+l)"  exp^-y)  ,94) 

where 

y  =  7(2  b+1 )/ r9b  2(b+l ) 2 1 

After  an  appropriate  value  of  the  lessen  ratio  )f  2  t: 
chosen,  one  iterates  (92)  ana  ^94)  to  fina  b.  The  ratio  oust  ue  :nosen 
empirically,  on  the  basis  of  net  ionizat’on  ana  -econo -nat 'on  -ates. 
value  of  this  ratio  which  yields  good  answers  is  0.1. 

The  next  step  is  to  find  the  normalized  population  distriou- 
tion.  A  Boltzman  distribution  is  a  gooa  approx -mat  ion  ‘or  tnese  .tites 
because  the  rate  of  col  1 isional  diffusion  of  electrons  is  very  hign  among 
very  highly  excited  states  compared  to  the  rate  of  transition  from  such 
states  to  lower  ones.  So  we  expect  tne  individual  state  ono.'r  r  *  - 
some  state  k  to  be 


r.  -  •  /.u. 
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=  C'g^e  *  =  C  k"  exp  rV/\  ek‘)  '  ,95, 

for  ic  between  D  and  the  maximum  (from  3.3.3).  to  find  tne  normali¬ 
zation  constant,  C,  eq.  (95)  must  be  summed  and  set  equal  to  tne  total 
population  of  the  artificial  state,  p  .  As  might  be  expected  from  the 

<3 

form  of  (95)  the  Summation,  even  in  an  integral  approximation,  is  pretty 
messy.  An  approximate  function,  accurate  to  better  than  4%  everywhere, 

C  =  2  p  (  e/'V )  3/  2/F  ,96: 

d 

with 


x  (b }  , ,  .  ,  , 
•  e  -  '  S  ( x  i  o ) ) 


S  ■.  x  (  $ 


i  i  /  e 


Radiation  Oscillator  Strengtn.  To  find  tne  effective  oscillator 


strength  *j r  relation  to  the  artificial  state  from  sane  real  state  t  we 
n j s t  now  sum  .46.'  over  all  states  which  make  up  the  artificial  state. 
Substitution  of  (39),  v90),  v 95 } ,  and  (96)  into  (48)  and  summing  yields 


!  =  "  r  (k,0  »  4A  V  2(  e/v )  3/  2  p  /(IF)  r  llP-LV/ik-J 

r  K  Sp  d  P  k(  k  2-  t2) 


(100) 


If  the  son  in  (100)  is  replaced  by  an  integral  its  value  is 


I  ;  ”  exp(  V/  9k  2)/  ’k(k  2-  i2)  dk 
D 


=  exo  v  *i\-  /  m  l2)  rE,  'i  (i_  -  1.)  '  -  E,  -  (i-  -  i_) 


(101) 


/•'■e-e  I .  s  tne  exponential  ’ntegrji  )f  tne  nrst  <inu.  It  is  convenient 
"nit  1  ,0  .'an  oe  approx ’mated  analytical?/  out  since  the  exponential 
ntegra’s  wi ' 1  oe  approx ’mated ,  there  is  ample  opportunity  for  truncation 


err.jr,  *o  nand’e  this  Je‘"’ne 


h  5  expi  -V/  -*o‘)  I 


■  A 
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V.-. 


1 108) 


-’■A*. 

m 


rhen  evaluate  h  for  two  cases.  Define 


y  ;  V/( “H; 


(103) 


•\  A 


i-  a 


H  =  exp(  a)/(2n2)  [tng/a  +  yU^  +  a2(cd-e)x2  +  a3(  a2+ct6+B2)x3 


+  a4y(  o^s)  2X4] 


(104) 


where 


at  =  .99999193,  a2  =  -0.24991055,  a3  ■  .05519968,  a4  =  -.00976004 
as  l/l2  -  1/ b 2  ,  8  =  l/l2  -  1/S2  ,  y  =  1/S2  -  1/b 2 


Case  2  y  >  1 

If  i2/ b2  <  10~2  Ex(y)/y  then 

H  =  i.  (1/b2  -  1/S2)  exp(-V/9D2) 


(105) 


If  i 2 / b 2  ^  10"2  £x(y)/y  then 

H  =  (  9/2  V  t)(  rEx(Va/9)/a  -  exp(-Vy/9)  Ex(  V  S/  9)  /  3  ]  (106) 

where,  from  Ref.  5, 

Ex(x)  =  (a2+a1x+x2)/(b2+b1x+x2)  (107) 


=  2.334733,  a2  =  0.250621,  b^  3.330657,  b2  -  1.681534 

Then  substitute  (102),  and  (101)  into  (100)  to  find  R.  Then 
equate  R  to  (48)  evaluated  for  states  i  and  a  to  obtain  at  last,  for  raai 


fR  =  2 v1/2e3/2  g.  exp(V/9b2)  H/  f(E  -E  )  213F  ] 


(108) 


where  g^  is  taken  from  (7) 

5. 7.2.2  Electron  Collision  Oscillator  Strength.  To  find  the  effective 
oscillator  strength  for  electron  collisions  to  the  artificial  state  from 
real  state  i,  sum  (52)  over  all  states  which  make  up  the  artificial 
state.  Substitution  of  (89),  (90),  (95),  and  (96)  into  (52)  and  summing 
yields 

Sm 

Re  =  Co  l3  ^  1(5  exP [V/(ek2)  ]/(k2-a2)4 

b 

Cq  i  2CECNe9Pagi/rV5/2Fl  (109) 

The  sum  in  (109)  is  also  rather  messy  and  gives  an  unrealis¬ 
tically  large  rate  when  b  approaches  i  (that  is,  when  9  and  are 
almost  out  not  quite  high  enough  to  drive  the  artificial  state  into  equi¬ 
librium  with  the  uppermost  real  state).  This  is  not  shocking  since  the 
assumption  of  zero  populations  for  states  below  b  and  full  Boltzmann  popu 
lation  above  is  over  simple.  An  approach  which  largely  alleviates  this 
problem  is  to  neglect  contributions  from  the  first  two  terms  in  the  sum; 
repl ace  b  with 


B  =  b  +  2 


(110) 


The  sum  in  (109)  can  be  approximated  by  replacing  it  with 


Replacement  of  k  by  ( V/ e) (1/ t2-l/k 2)  as  the  variable  of  integration 
reduces  the  integral  to  a  fourth  order  exponential  integral,  which  can  be 
reduced  to  a  first  order  exponential  integral  by  integration  by  parts.  T 
avoid  truncation  problems  the  following  cases  are  recognized. 

Case  1  V/(eB2)  <_  0.1  and  4z2/B2  £  0.1 

Then  an  adequate  approximation  is 

R  =  (l/b2-l/S 2) /2  (112) 

Case  2  Case  1  failed  but  (V/e( (l/e2-l/B2)  <  10.  Define 
8  i  (V/e)  (1/ 12  -  1/B2) ,  T  i  (v/e)(l/t2-l/S2) 

R  =  ( 1 Z  e 2 )  ~  1  rV/(-u2}13  -'exp  r(V/  oB2)  ’( 1-1/ *+2/ j2  -  I  x  (  0  )  /  } 

-  exp  'V/(  eS2)  .1  -  I /  r  2/  r:-Ex^  t ) )/  rl  Hi; 

where  tne  function  Ex:xl  is  that  of  ,107). 

Case  3  Not  Case  1  )r  Case  2 

R  3  (2  t2 )"  1  rv/(  9 1 2 )  13  lex  p  r ( V / ^B  2 )  3  (1  -  4/  «  ♦  20/ b  2 
-  1 20/  ? 3  +  540/u14  -  6720/ ♦  60,480 '9*0 


SECTION  6 
BASIC  DATA 


In  addition  to  normal  mathematical  and  physical  constants,  the 
major  data  requirements  of  this  calculation  are  quantities  associated  with 
energy  levels  and  especially  oscillator  strengths  for  transitions  among 
energy  levels. 

Energy  levels  and  their  associated  principal,  orbital,  and  spin 
quantum  numbers  are  obtained  primarily  from  Ref.  7.  Oscillator  strengths 
are  obtained  primarily  from  Ref.  8. 

A  significant  amount  of  art  's  involved  in  this  process,  while 
not  all  needed  data  exists,  particularly  for  oscillator  strength,  far  more 
energy  ’eve's  are  vnown  than  :an  oe  nooeiled.  It  is  necessary  to  comoine 
n any  states  togetner  to  reduce  the  total  number  noueled  to  scout  10.  In 
io  i  ng  so,  one  tries  to  combine  states  whose  energy  1  eve’s  are  nearby, 
avoid  creating  jr 1 1 f i c i a  1  japs  in  the  distribution  of  energy  levels,  not 
eve's  with  different  ,ptns,  nor  j> *•  >'eront  pr  <nc  ’  pa!  duanf  jn  njmber. 
Occasionally  conflicts  arise  among  *;nese  criteria  and  are  ••esolveu  by  best 
jjdgenent.  nnotner  proolem  ar'ses  wnen  ail  egenstates  if  i  given  pr -nc 
pal  quantum  number  are  not  described  in  the  literature.  In  that  case  we 
jse  available  data  tc  letermine  node’ed  state  energy,  <>tc.,  and  r^ise  *he 
statistical  weight  to  account  for  the  missing  states. 


Oscillator  strengths  are  summed  over  upper  levels  and  averaged 
over  lower  leve's.  For  transitions  allowed  by  LS  coupling  rules  but  which 
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In  Griem’s  formula  for  Stark  broadening  (Ref.  2),  the  energy 
gaps  to  the  next  lower  and  next  higher  orbital  state  appears.  Since  we 
frequently  lump  states  of  differing  orbital  angular  momentum  together, 
these  gaps  are  inserted  as  data,  averaged  as  best  as  we  are  able. 

An  example  of  the  results  of  this  process  is  given  in  Table  1 
which  shows  data  entered  into  the  computer  for  singly  ionized  oxygen.  In 
this  case  63  eigenstates  from  Ref.  7  were  reduced  to  10.  For  the  upper  5 
model  states,  preservation  of  spin  was  given  up,  as  can  be  seen  from  the 
non-half  integer  values. 
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Table  1.  01+  Data.  The  first  column  is  state  index,  k,  the  second  state 

energy  (cm-1),  then  designation,  statistical  weight,  g,  princi¬ 
pal  quantum  number,  n,  orbital  quantum  number,  i ,  spin  quantum 
number,  s,  mean  energy  difference  to  i-l,  w-  (cm-1),  energy  dif¬ 
ference  to  z+1,  w+  (cm-1),  and  oscillator  strengths  where  the 
number  in  parentheses  is  the  value  of  k  for  the  upper  state. 
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SECTION  7 
SAMPLE  RESULTS 


To  date  little  useful  data  for  direct  comparison  of  overall 
results  has  been  found.  Necessary  basic  comparisons  have  been  made.  For 
example,  our  fit  to  Lotz's  curve  has  been  checked;  Lotz  has  compared  his 
curve  to  electron  collisional  cross  section  data.  And  we  have  checked  our 


photon  ionization  cross  sections  against  Kramer's.  This  sort  of  fundamen¬ 
tal  comparison  is  required  but  accomplishes  no  more  than  assurance  that 
the  formulae  were  correctly  encoded;  such  checks  do  not  give  evidence  of 
the  adequacy  of  the  method  to  produce  correct  answers.  The  only  way  to 
check  the  total  method  is  to  compare  to  other  calculations  and  to  limiting 
cases  over  as  wide  a  range  as  possible. 

7.1  COLLISIONAL-RADIATIVE  RECOMBINATION  AND  IONIZATION. 


The  most  severe  test  of  the  method  we  have  found  is  based  on  the 
fact  that  a  complete  temporal  calculation  of  electron  eigenstate  popula¬ 
tions  lust  automatically  account  for  ooth  radiative  ano  collisional  ioni¬ 
zation/recombination  processes.  The  only  prior  calculations  known  which 
are  in  a  form  for  ready  comparison  are  those  in  Ref.  9  by  Sates,  Kingston, 
and  McWhirter  (hereafter,  BKM) .  The  BKM  calculations  are  restricted  to 
hydrogen  or,  separately,  nydrogenic  ions  in  an  optically  thin  nedium*, 
using  an  electron  collisional  cross-section  derived  from  Grysinski's  semi 
classical  ionization  rate  coefficient  to  back  out  state- to- state  excita¬ 
tion  cross  sections  for  electron  collisions. 

*  They  also  carry  out  a  version  of  a  calculation  for  a  medium  optical 1  y 
tnicx  to  series  of  i ines.  They  jo  not  seem  to  consider  tie  ’on  ’  z'iu 
effect  of  the  radiation,  we  have  not  yet  attempted  to  compare  to  tn's 

work . 


Their  calculations  make  the  following  major  assumptions,  which 
in  the  following  will  be  referred  to  by  the  letters  in  parentheses: 

(P)  Only  the  ground  state  is  significantly  populated 

(D)  Only  the  ground  state  requires  a  differential  equation  solution, 
all  others  can  be  approximated  as  in  quasi-steady  equilibrium. 

(E)  All  energy  levels  are  present  at  all  electron  densities,  that 
is,  no  effect  of  plasma  polarization  is  included. 

Over  the  range  of  applicability  of  these  assumptions,  BKM  expect 
their  results  to  be  accurate  to  about  a  factor  of  3.  Our  method  is  more 
general  in  that  it  makes  none  of  the  three  assumptions,  P,  D,  or  E  but 
less  accurate  in  that  for  practical  reasons  we  are  limited  to  principal 
quantum  number  of  3  or  4.  To  check  our  method  we  ran  hydrogen  using  only 
the  first  four  eigenstates  as  well  as  the  first  10  eigenstates.  If  our 
nethod  were  perfect,  we  would  get  the  same  answers  in  both  cases. 

7.1.1  Recombination  Coefficient  for  Hydrogen. 

Figures  1  througn  5  compare  our  calculated  values  of  col  1 isional 
'aoiative  -ec omo mat  ion  coefficient,  a,  with  that  of  BKM.  In  all  figures 
a  vertical  Qar  indicates  the  value  of  electron  density,  N  ,  at  which 
one  or  more  BKM  assumptions  Degin  to  fail,  for  higher  the  SKM 
approx imation  is  not  reliable.  In  Figures  1,  2,  and  3  assumptions  D  and  E 
fail  first,  in  Figures  5  ana  6  assumption  3  fails  first. 


It  can  oe  seen  that  wnile  our  method  is  not  perfect,  the  4  state 
case  is  pretty  consistent  with  the  10  state  case  and  both  are  consistent 
w’tn  BKM  over  the  range  of  for  which  the  BKM  approx imation  holds. 

In  *ne  -  ad!  at  we  recombination  ’imit,  at  low  N  ,  our  a  is  consistently 
below  BKM  oecause  the  i_otz  rate  is  oelow  that  jsea  oy  BKM,  Dut  is  at  least 
dose  tn  the  factor  of  3  BKM  expects. 


sVV*v*r 


Figure  4.  Hydrogen  Comsional  Radiative  Recombination  Coefficient  at 
16,000*K.  The  8KM  coefficient  is  shown  as  a  solid  line.  Our 
values  for  10  states  are  shown  as  open  circles,  for  4  states  as 
open  squares.  BKM  is  not  reliable  to  the  right  of  the  vertical 
line. 


■  -  v. ^■.V.'^'AV/aVA-.VoV^.  -\V T A.  A  < -  V,' <.V.V .V .'  r 


S  's*  4 


«C 


r 


’34,jOC**.  J**‘  <?«*  r  *, 

*a\jes  for  ij  states  are  snowr  as  ,per  r, 
ope*i  squares,  s**  ’s  o(,t  r*  af-  e  * .  *ie  «- 
'  1n<  rX/  ’  Q  '. 'r«.  ’  es  ano  squares  ,s#-  {**>  je 


At  tne  two  lowest  temperatures ,  Figure  1  and  2,  the  4  state  cal- 
c/at'an  setters  a  glitch  at  about  the  value  of  for  which  BKM 
degrades.  This  is  almost  certainly  a  problem  in  our  effective  electron 
collisional  oscillator  strength  from  low  lying  states  to  the  artificial 
state,  due  to  tne  fact  that  the  gap  in  energy  between  state  4  and  the 
artificial  state  is  aoout  40  times  the  electron  temperature.  This  gap  is 
'-educed  almost  a  factor  of  6  for  the  10  state  case,  and  of  course,  is  also 
reduced  at  higher  electron  temperature.  We  do  not  consider  this  a  serious 
yrob’em  because  we  ao  not  expect  to  encounter  such  high  at  such  low 
e'ectron  temperature.  If  tnis  expectation  is  wrong,  it  will  be  necessary 
ts  either  improve  the  method  or  carry  more  states. 

At  high  electron  temperature,  especially  Figure  5,  our  a  climbs 
it,’..*  tnat  ;f  BKM  at  high  .  In  part  this  is  due  to  a  differ- 
e*  c *  'r  jur  a e r  i n  1 1 1 o n  of  a  ana  tneirs.  We  count  recombination  to  any 
!'.  ,.<v.  ,-r  t,  wnereas  3  KM  count  only  recomD  ination  to  the  ground 
*n*r’  t  -M  •  ••  js>,.Mipt  ion,  is  /alia  the  two  definitions  are  numeri- 
:  .  j  e"t ,  out  wnen  the  upper  states  acquire  a  significant  popul  a- 

■  ,  .  ,  v  .e  ‘  j '■  tne  rase  )f  cigure  5  at  nign  they  are  very 

.  •• . *.  «».•  t  o ‘nation  •;  aef  f  ic  lent  to  the  ground  state,  to  maxe  a 

-  -■  ^  m,  s  jfiown  -n  -igure  5  as  so!  id  circles  ana 

.  r  >-  m  -nat  the  agreement  is  quite  good,  indicating  tnat 
fw  a  /a’ ion  is  less  sensitive  to  assumption  ?  than  mignt  De  imag- 

•  r 

'.1.2  Ion i /at ion  Coefficient  for  Hydrogen. 

;  j, n,  -  and  h  snow  fxnpar i sons  of  our  co  1  1  i s i onal - r ad ’ a t i  ve 
■  /a’  '  *♦*''  ' **nt ,  . ,  to  those  of  BKM.  Our  values  are  uncomfortably 
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Figure  6.  Hydrogen  Collisional  Radiative  Ionization  Coefficient  at 

4000’K.  The  8KM  coefficient  is  shown  as  a  solid  line.  Our 
values  for  10  states  are  shown  as  open  circles,  for  4  states  as 
open  squares.  BKM  is  not  reliable  to  the  right  of  the  vertical 
1  ine. 
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Figure  7.  Hydrogen  Col  1 1  $ lonal  Radiative  Ionization  Coefficient  at 

16,0Q0*K.  The  BKM  coefficient  is  shown  as  a  solid  line.  Our 
values  for  10  states  are  shown  as  open  circles,  for  4  states  as 
open  squares.  8KM  is  not  reliable  to  the  right  of  the  vertical 


Figure  8.  Hydrogen  Col  1 islonal  Radiative  Ionization  Coefficient  at 

64,000‘K.  The  BKM  coefficient  is  shown  as  a  solid  line.  Our 
values  for  10  states  are  show  as  open  circles,  for  4  states  as 
open  squares.  BKM  is  not  reliable  to  the  right  of  the  vertical 
1  ine. 


i '  a'  it:  ins  nidifying  several  aspects  of  the  artificial  state  show 
•'M1  *  ne  »  i-  4e  it  ,  at  4uu0  V  is  sensitive  to  such  details.  For  example, 
if  tne  pe’arization  depression  energy  of  £q.  [3)  were  increased  a  factor 
jf  4,  essent'a’.ly  all  o  i  screpanc  ies  in  S  left  of  the  vertical  line  would 
disappear  for  all  three  figures,  6-8.  Although  this  would  be  the  simplest 
way  to  eliminate  tne  discrepancies,  we  don’t  feel  qualified  to  argue  the 
accuracy  of  v 3 ) ,  besides,  there  are  other  options  for  eliminating  the  dis 
crepancy.  Because  S  is  so  small  at  400Q’K,  we  don’t  feel  it  worth  the 
effort  to  complicate  the  code. 

More  serious  discrepanc ies  occur  in  the  region  where  BKM  breaks 
down.  Our  value  of  S  tends  to  be  higher  than  BKM  in  this  region  because 
( 1 )  polarization  depression  of  ionization  energy  makes  ionization  more 
favoraole  when  electron  temperature  is  less  than  ionization  potential  and 
as  population  of  upper  states  oecomes  important,  either  in  steady 
state  jr  transiently  luring  recomo mat  ion ,  ionization  of  these  states  :an 
lom'nate.  'he  'irst  if  tnese  effects  occurs  wnen  BKM  assumption  £ j  s 
/•plated,  tne  second  can  occur  transiently  wnen  them  assumption  /J  l  *s 
violated  anq  permanently  when  them  assumption  •;?)  is  violated. 

'ne  scatter  ;n  jur  "'ates  at  very  nign  occurs  Jue  to 
breakdown  of  assumption  ( D ) .  when  the  quasi-steady  assumption  fails  for 
several  excited  states  over  a  significant  fraction  of  the  time  necessary 
for  the  process  to  go  to  completion,  neither  a  no r  S  is  a  constant.  The 
particular  value  ootained  depends  on  the  particular  moment  chosen  for 
evaluation.  The  range  of  values  obtainable  is  great,  as  illustrated  by 
the  orders  of  magnitude  difference  between  S  at  Ng  =  102b  from  Figure 
7  and  Figure  8,  and  by  the  difference  in  4  state  and  10  state  values  in 
Figure  8  at  Ne  =  1026. 

In  summary,  where  the  value  of  S  is  significant  ano  the  BKM 
assumptions  valid,  agreement  is  sati sf actory .  Where  S  is  uninterestingly 


CO*l- 


sma’I,  our  values  are  nigh,  where  BkM  assumptions  tail,  the  normal 
sept  of  S  tireass  down  and  no  calculations  are  available  for  c  oinp  ar  i  son . 

1.2  APPROACH  TO  STEADY  STATE  WITHOUT  RADIATION. 

In  addition  to  comparison  to  the  results  of  BKM,  one  can  con¬ 
sider  the  "reasonableness"  of  the  final  steady  state  configuration  for  non 
equilibrium  conditions  and  the  rate  and  manner  of  approach  to  steady 
state . 


A  fairly  serious  problem  with  any  such  discussion  is  finding  a 
str aightforward  and  clear  method  of  presentation  of  results  for  a  system 
as  complex  as  several  dozen  eigenstates  distributed  over  several  states  of 
ionization  for  an  element  of  interest  such  as  nitrogen  or  oxygen. 

In  an  attempt  to  invent  an  appropriate  graphic  technique,  we 
nave  adopted  the  following. 

Consider  the  partition  function  of  some  state,  <,  of  an  element 

at  temperature  *  and  electron  density  N  . 

J  e 

p(k)  =  (Qe/Na)Z  gk  expt-Ek/e)  ,116) 

where  Qe  is  the  partition  function  for  free  electrons, 

0e  =  6. 03 7x10 21  93/:  (cm*3)  (117) 

e(eV)  is  temperature,  Ng  (cm-3)  is  electron  density,  gk  is  statistical 
weight  of  the  state  k,  is  energy  of  state  k  relative  to  the  lowest 

eigenstate  of  the  element  (normally  the  ground  state  of  the  neutral,  but, 
e.g.,  for  oxygen,  the  ground  state  of  J" ) ,  and  Z  is  the  charge  state  of 
the  ion  to  which  eigenstate  k  belongs. 


ne  total  partition  function  of  the  element  is 


In  equilibrium,  the  fraction  of  the  element  in  state  k  is  given 


f ( k )  =  p(k)/P  =  r(Qe/Ne)  9k/p  1  ex p ( -E k / 9) 


(119) 


F(k,9,Ne)  f ( k )  =  exp(-Ek/9) 


(120) 


F(k,  Me:  =  =/  riQe'Ne)^  gk 


The  form  (120)  is  a  rather  versatile  way  to  present  results. 

One  can  graph  the  exponential  of  the  RHS  against  E^  to  obtain  a 
straignt  line  on  semi-log  paper,  calculate  F  a  priori,  then  use  whatever 
f( k)  results  from  a  calculation  to  evaluate  the  LHS  and  graph  it.  If  the 
calculated  f(k)  are  in  equilibrium  at  9,  the  two  curves  will  lie  on  top  of 
each  other,  if  not,  one  has  a  visual  illustration  of  the  degree  of  depar¬ 
ture  from  equilibrium.  The  disadvantage  of  the  presentation  is  that  a 
curve  of  the  LHS  of  (120)  is  effectively  divorced  from  the  actual  values 
of  f(k),  so  that  no  clue  exists  as  to  which  electron  or  ionization  states 
are  most  populous. 
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7.2.1  Time  Development  of  Hydrogen. 

An  introduction  to  this  sort  of  presentation  is  provided  by  Fig¬ 
ure  9  which  shows  the  results  for  hydrogen  initially  neutral,  bathed  in 
108  electrons  at  1  eV  with  no  radiation  field.  Our  initialization  over¬ 
populates  the  ground  state  (k=l)  and  the  first  excited  state  (k=2);  as 
shown  by  the  aashed  line,  relative  to  the  1  eV  solid  line.  In  slightly 
more  than  104  sec  the  states  have  come  to  a  steady  configuration  with 
electron  collisional  processes  biased  by  spontaneous  radiation,  shown  by 
the  dotted  line.  In  this  case  one  might  speak  of  an  "effective  eigenstate 
temperature"  of  about  0.4  eV  to  approximate  the  electron  distribution.  In 
the  final  state  the  actual  fractional  populations  are  0.945  for  state  1, 
6.6xl0~6  for  the  artificial  state  at  about  13.6  eV,  and  0.055  for  H1+ 
at  13.6  eV.  All  tne  others  range  from  10" 12  to  10~13. 

7.2.2  Hydrogen  Steady  State  vs  Electron  Density. 

Figure  10  shows  the  final  configurations  for  hydrogen  in  a 
plasma  of  1  eV  electrons  at  several  values  of  electron  density.  By  Ne 
=  1018  electron  collisional  processes  dominate  spontaneous  radiation 
sufficiently  that  an  equilibrium  distribution  is  attained. 

7.2.3  Time  Development  of  Nitrogen. 

A  much  more  interesting  and  much  more  complex  case  is  illus¬ 
trated  in  Figure  11.  Two  calculations  are  illustrated.  Both  have  elec¬ 
tron  density  of  1014  cm" 3  and  temperature  3  eV.  The  first  was  initial¬ 
ized  with  all  the  nitrogen  stripped  to  7+,  and  the  second  was  initialized 
with  all  nitrogen  neutral.  By  2.67xl0~4  sec  the  first  case  has  recom¬ 
bined  enough  to  extend  an  erratic  pattern  of  electron  population  down  to 
N4+,  by  1.23x10" 3  sec  populations  extend  to  2+.  By  a  few  tenths  of  a 


F(k »  6.N  ) f (k ) 


Nitrogen  in  10  *  electrons  at  3  ev.  The  heavy  solid  line  is 
the  final  steady  state.  The  lignt  uroken  line  extending  down 
to  4>  shows  states  at  2.67x10"“  sec  after  initialization  with 
only  7+  nitrogen  present.  The  light  solid  line  extending  down 
to  2+  shows  states  after  1 . 23x10" 3  sec.  The  dashed  line 
extending  through  1+  shows  states  1.87x10"  11  sec  after  initial¬ 
ization  with  only  0+.  Both  cases  approach  the  heavy  solid  line 
with  total  ion  populations  indicated  by  horizontal  lines  over 
the  energy  range  corresponding  to  the  state  of  ionization. 


•  -  •  .  •  <■-  »:  tif  so  >10  me,  which  is  better 

i  ■  •  >•-  -*  J  e»  t n an  tne  true  electron  tem- 

'-r  -e  ase  snows  popu'ation  through  1+  by 

.  )■  j  ,  --a.nes  iteaov  state  at  tne  heavy  solid  line  by 

’■m  -e  at'.e  pop./at’on  of  states  show  much  more  structure  than 
.  •  ?->.  a„se  oscT  ’ator  strengths  connecting  the  nitrogen  states 

.  ••  .  ■■  oon-'-ast  to  the  smooth  values  appropriate  for  fully 

ire  -lyjr^er  states.  Note  that  the  first  excited  state  of  neutral 

-  •  ue-  '  >  moe’-popu'  atea  by  nearly  4  orders  of  magnitude  relative  to  the 
>e  st  ate  at  I.S"*iO*1'  sec  for  the  case  initialized  with  neutral 

■  ■ - . ue” .  ’n>s  -s  Jue  to  the  fact  that  this  state  is  optically  forbidden 

-  •  -  • •  *'  * ne  moma  state.  Our  algorithm  arbitrarily  assigns  an 
e**ect'*r  osci'’ator  strength  of  10'3  to  such  transitions  for  electron 

-se  t  wouij  oe  even  more  severely  underpopulated. 

'ne  effective  ionization  rate  coefficient  from  N0+  to  N1+  for 
t-rs  case  is  aDOut  7  <10"  10  cm3/sec,  about  twice  the  BKM  value  for  hydro- 
:er-c  'ons  ana  about  seven  times  Sappenf i e 1 d 1 s  calculation  for  ground 
state  'onization  coefficient.  The  agreement  with  BKM  is  quite  satis¬ 
factory,  particularly  -in  view  of  the  fact  that  nitrogen  is  not  especially 
hydrogenic.  The  discrepancy  with  Sappenf ield  is  presumably  due  partly  to 
differences  in  collision  cross  section  and  also  is  due  in  part  to  contri¬ 
butions  from  excited  states,  especially  states  2  and  3,  which  acquire 
fairly  large  populations  (0.6  of  state  1  and  0.2  of  state  1  respec¬ 
tively). 

7.2.4  Nitrogen  Steady  State  vs  Electron  Density. 

Figure  12  shows  steady  state  configurations  for  nitrogen  in  a 
bath  of  1  eV  electrons  and  in  the  absence  of  a  radiation  field,  the  opti¬ 
cally  thin  limit.  When  electron  density  gets  as  high  as  1016/cm3 


Figure  12.  Nitrogen  at  1-eV.  Nitrogen  steady  state  eigenstate  distribu¬ 
tion  for  1-eV  electrons  and  no  radiation  is  shown  for  three 
values  of  electron  density.  Solid  circles  correspond  to  1018 
electron/cm3  or  higher,  open  triangles  to  10 12  electrons/cm3 
and  open  circles  to  108  electrons/cm3. 
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eigenstate  populations  follow  the  equilibrium  exponential  out  to  20  eV  or 
30  eV  before  falling  below  it.  By  electron  density  1018  cm" 3  all  states 
retained  by  the  code  are  within  1%  of  equilibrium;  states  above  1+  are 
deleted  as  described  in  Section  3.3.4  at  this  high  electron  density. 


Steady  state  distribution  of  nitrogen  eigenstates  for  electrons 
at  10  eV  are  shown  in  Figure  13.  Note  that  the  abscissa  is  logarithmic, 
so  the  range  goes  from  10'7°  to  109°.  The  curves  show  the  extreme 
sensitivity  to  electron  density  for  weak  plasmas,  in  that  twelve  orders  of 
magnitude  in  electron  density  (10~2  to  1010)  result  in  12  to  40  orders 
of  magnitude  difference  in  state  population  relative  to  equilibrium.  It 
also  illustrates  the  disadvantage  of  this  presentation;  the  ground  state 
of  N0+  is  overpopulated  by  an  impressive  factor  of  1093  but  even  so,  all 


states  account  for  only  1  part  in  10 12  of  the  nitrogen. 


The  approach  to  equilibrium  is  shown  in  more  detail  in  Figure 
14.  Here  it  can  be  seen  that  even  1018  cm" 3  leaves-  the  lower  states 
about  an  order  of  magnitude  overpopulated  wnereas  102°  cm" 3  brings  them 
into  equilibrium.  If  the  curves  of  Figure  14  were  to  be  extended  to 
higher  energy,  they  would  show  that  even  at  102°  cm"3,  the  excited 


states  of  N3  around  700  eV  to  800  eV  would  be  severely  underpopulated. 


7.2.5 


Value  of  Electron  Density  for  Equilibriun  in 
Optically  Thin  Plasma. 


For  electron  collisions  to  drive  the  eigenstates  into  equilib¬ 
rium  in  an  optically  thin  plasma,  it  is  clear  that  electron  collisional 
excitation  must  dominate  spontaneous  radiation.  Given  a  ladder  of  some 
dozens  of  eigenstates  one  would  expect  the  downward  bias  to  accumulate,  so 
electron  collisions  might  need  to  dominate  spontaneous  radiation  by  one  or 
more  orders  of  magnitude. 


EIGENSTATE 


The  ratio  of  electron  collisional  excitation  from  lower  state  i 
to  upper  state  k  to  spontaneous  radiation  from  k  to  i  is,  from  (48)  and 
(51) 


Rp  =  Ne^gk/9t^Cec/Ar^  exp  [-(Ek-E  ?)/e]/e1/2(Ek-E  f) 3 


1.5x10- 13  exp [-(E.-E  )/ e]/ [e1/2(E.-E  ) 3 ]N 


(122) 


when  e  and  E  are  measured  in  eV  and  N  in  cm*3,  and  the  ratio  of 

e 

statistical  weights  has  been  taken  to  be  unity.  Then  a  distribution  ap¬ 
proaching  equilibrium  can  be  anticipated  whenever  Rg  »  1,  or 


N  >  6xl012el/2  (Ek-E  ) 3  exp[(E,-E  )/ e] 


(123) 


If  Ek~Ez  -  10  eV  as  appropriate  for  lower  ionization  states  of 
nitrogen,  or  for  neutral  hydrogen,  one  needs 

>  5  xlO 16  cm  3  . 


Thus  one  should  not  be  surprised  to  see  hydrogen  states  come  into  equi¬ 
librium  at  about  1018  cm* 3  and  the  nitrogen  line  of  Figure  14  parallel 
to  the  equilibrium  line  for  the  first  few  states  of  ionization  at  1018 


On  the  other  hand,  if  E.-E  -  100  eV  as  is  the  case  for  some 

k  i 


states  of  N4+  and  most  of  N5+,  one  needs 
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IV  >  5x10  J  cm 


for  equality.  These  states  aren't  likely  to  be  brought  into  equilibrium 
for  Ng  less  than  about  1025  cm*3.  Long  before  this  electron  density 
is  reached:  (1)  equilibrium  populations  of  highly  ionized  states  are 
negligible  (<■  10_2°  for  N4+)  and  (2)  our  simple  treatment  of  polariza¬ 
tion  effects  has  become  inappropriate. 


7.3 


APPROACH  TO  STEADY  STATE  WITH  PLANCK  RADIATION. 


The  rate  of  excitation  for  an  allowed  bound-bound  transition  of 
energy,  e,  by  a  Planck  spectrum  is 

rBB  =  f  a(  v)  U(v)dv  (124) 

If  the  line  is  reasonably  narrow,  as  is  almost  always  the  case,  the  Planck 
spectral  flux,  U(v),  can  be  removed  from  the  integral  in  (124)  and  the 
integral  evaluated  (Ref.  3,  p.  59), 

oo 

rBB  “  ^(v)  /  tf(u)dv  =  f  U(  v)  (125) 

—  oo 

where  rg  is  the  classical  rauius  of  the  electron,  c  the  speed  of 
light,  and  f  the  oscillator  strength. 

From  (60),  'J(v)  can  be  written  as 

U(v)  =  Cp€2/ (exp(  d  9)-l )  (126) 

Then  for  a  transition  of  energy,  e,  and  nominal  oscillator 
strength  of  0.1,  the  time  to  reach  equilibrium  is  of  order 

T  =  r-1  =  2*10- 7  (exp(e/9)-l)/e2  (127) 

BB  BB 

Thus,  for  a  Planck  spectrum  at  1  eV  temperature,  we  should 
expect  equilibrium  among  bound  states  in  about  10"5  sec  if  the  biggest  gap 
between  states  is  about  8  eV,  which  is  about  the  case  for  N0+,  N1+,  and 
N2+,  or  about  1.5  sec  for  a  biggest  gap  of  22  eV  (N3+). 


To  bring  the  states  of  ionization  into  equilibrium  is  more  dif¬ 
ficult.  The  cross  section  for  photoionization  is  given  by  (67).  If  this 
value  is  substituted  into  (124)  with  proper  care  for  units,  one  finds 

rBF  =  CHCpIZ/n  /  Mexp(e/e)-l)]-1  d£ 

h 

-  CHCpIze/n  exp(-I2/e)  (128) 

or 

tbf  =  rg,J  -  2x10“ 7  exp(I2/e)/(Iz0)  (129) 

where  (129)  has  assumed  that  Iz  >  0  and  n  =  3. 

If  one  assumes  bound  states  to  come  into  equilibrium  oefore 
ionization  states  do,  then  we  can  apply  (129)  to  the  ground  state  and  be 
within  a  few  orders  of  magnitude  of  the  right  answer,  since  excited  states 
lose  in  population  all  that  they  gain  by  the  Boltzmann  factor  in  (128). 

Application  of  (129)  yields  3"2  sec  for  N0+  to  N1+,  5xl04  sec 
for  1+  to  2*,  2x10 12  sec  for  2+  to  3+,  and  1025  sec  for  3+  to  4+  for  a 
1-eV  Planck  spectrum.  These  numbers  verify  the  assumption  made  above  that 
eigenstates  within  a  given  ion  should  equilibrate  faster  than  different 
ions. 


Figure  15  shows  nitrogen,  initially  neutral,  bathed  by  a  1.0  eV 
Planck  spectrum  at  2*10* 5  sec,  25  sec,  and  IQ21  sec.  To  reduce  effects 
of  electrons  to  a  minimum,  electron  density  was  taken  to  be  10" 10 
cm"3,  and  to  ensure  that  any  appreciable  collisional  effect  would  show, 
at  least  on  the  steady  state  curve,  electron  temperature  was  set  to  0.1 
eV.  One  can  see  that  eigenstates  of  0+  and  1+  (the  only  ions  present'  are 
pretty  well  in  equilibrium  by  2x10" 5  sec,  as  expected  from  (127)  for 
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25  sec 


bound-bound  transitions  (their  lines  are  parallel  to  the  1  eV  line).  By 
25  sec  the  0+  ion  state  and  1+  ion  state  are  almost  equil iberated  (lines 
are  parallel  to  the  1  eV  and  almost  continuous),  as  expected  from  (129). 

By  1021  sec,  all  states  are  equilibrated,  the  curve  is  indistinguishaole 
from  an  exponential  at  1  eV.  At  10 6  sec  (not  shown)  1+  and  2+  are  in 
equilibrium  as  expected,  while  2+  and  3+  are  still  far  from  equilibrium, 
as  expected.  The  fact  that  3+  equilibrates  with  4+  by  1021  sec,  rather 
than  requiring  1025  sec  is  a  bit  mysterious.  At  least  one  order  of  magni¬ 
tude  is  explicable  by  the  distribution  of  eigenstate  population  in  3+, 
which  makes  (129)  pessimistic. 

Figure  16  illustrates  a  situation  which  is  a  bit  more  interest¬ 
ing  physically,  electrons  hotter  than  radiation.  Nitrogen,  initially 
neutral,  is  bathed  in  a  1  eV  Planck  as  in  the  case  of  Figure  15,  but  the 
electrons  have  a  temperature  of  10  eV  and  densities  from  10'10  cm- 3  to 
1020cm-3.  Shown  are  final  steady  state  population  distributions.  For 
the  case  of  10“ 10  cm- 3  the  portion  of  the  curve  for  0+  is  an  artiface  of 
the  computer  plotting  program  and  should  t>e  ignored;  there  is  in  fact  so 
little  0+  present  (total  ion  fraction  less  than  10~23)  that  the  computer 
has  zeroed  out  all  0+  states.  The  negative  glitch  in  the  upper  states  of 
1+  should  likewise  be  ignored  -  one  state  has  too  little  population  to  be 
calculated.  The  other  irregul ari ties  in  the  10"10  curve  are  due  to  elec¬ 
tron  collisional  effects  on:  (a)  very  closely  spaced  states,  particularly 
that  from  the  artificial  state  to  the  ground  state  of  the  next  ion  and  (b) 
states  so  widely  separated  that  the  1  eV  Planck  spectrum  cannot  bridge  the 
gap  even  as  well  as  10”  10  electrons/ cm 3.  This  latter  is  the  case  above 
about  250  eV  eigenenergy.  As  electron  number  density  increases  the  curve 
continues  to  parallel  the  1  eV  radiation  temperature  over  some  ranges  of 
eigenstates  until  somewhere  between  1  electron  cm- 3  and  10 10  electrons 
cm-3  where  the  curve  shape  loses  definition.  Between  10  14  cm- 3  and  10 18 
cm' 3  the  electrons  take  over  and  by  102°  cm- 3  the  eigenstates  through  N4+ 
have  been  forced  into  equilibrium  at  10  eV,  despite  the  cooler  photon 
fl  ux. 


7.4 


STEADY  STATE  WITH  LINE  RADIATION. 


A  situation  of  considerably  more  physical  interest  than  the  pre¬ 
ceding  examples,  and  one  which  has  an  important  implication  for  the  accu¬ 
racy  of  previous  calculations,  is  that  in  which  line  radiation  is  present 
with  no  continuum.  To  investigate  the  effect  of  a  few  strong  lines,  we 
chose  an  electron  density  of  1012/cm3,  corresponding  to  high  enough  alti¬ 
tude  that  the  continuum  should  be  extremely  weak  and  not  too  many  strong 
lines  are  expected.  Electron  temperature  was  taken  to  be  10  eV,  a  typical 
value  expected  when  line  radiation  is  an  important  energy  loss  mechanism. 
Under  these  conditions  with  no  radiation  present,  initially  neutral  nitro¬ 
gen  approaches  steady  state  in  time  of  order  0.1  sec  with  state  4+  domi¬ 
nant;  fractional  populations  in  steady  state  are: 

0+  1+  2+  3+  4+  5  + 

6*10-12  5*10"7  2  *10" 3  0.2  0.71  .087 

Since  4+  is  the  dominant  ion  we  investigated  the  effect  of  add¬ 
ing  a  few  lines  which  connect  4+  states.  All  lines  were  taken  to  be  satu¬ 
rated  at  the  Planck  limit  for  10  eV.  They  were  added  in  order  of  increas¬ 
ing  energy. 

The  first  case  is  illustrated  in  Figure  17  which  shows  the 
,  cgion  of  eigenenergy  occupied  by  the  electron  states  of  N4+.  Two  curves 
are  shown,  that  for  no  radiation  and  that  for  a  single  line  of  2.69  eV  at 
the  Planck  limit  connecting  states  3  and  4  of  N4+.  The  cases  are  indis¬ 
tinguishable  except  that  the  line  has  raised  the  population  of  state  3 
about  a  factor  of  2.  This  is  not  an  impressive  effect,  but  perhaps  one 
should  not  expect  too  much  from  a  line  whose  intensity  is  over  an  order  of 
magnitude  down  from  the  peak  of  the  Planck  spectrum  and  which  connects 
states  whose  populations  are  rather  small  compared  to  the  ground  state. 


The  next  calculation,  shown  in  Figure  18  is  a  little  bit  more 
impressive.  Here  we  added  a  line  at  9.9978  eV  connecting  states  1  and  2. 
One  can  see  that  state  2  is  brought  up  to  approximately  the  correct  equi¬ 
librium  ratio  to  state  1  (the  line  connecting  state  1  to  2  in  Figure  18  is 
nearly  parallel  to  the  10  eV  line).  Further,  several  of  the  upper  states, 
which  have  optically  allowed  transitions  to  state  2,  have  increased  popu¬ 
lations.  This  is  the  sort  of  result  one  would  expect;  in  fact,  a  common 
method  of  assessing  the  effect  of  a  medium  optically  thick  to  the  line 
connecting  ground  state  to  first  excited  state  is  to  set  those  populations 

in  equilibrium  ratio,  then  solve  the  remainder  in  the  optically  thin 

approximation.  In  an  element  with  many  ionization  states  such  as  nitro¬ 
gen,  however,  the  line  has  another  effect  which  is  at  least  as  important 
as  redistributing  the  4+  state  populations.  A  line  with  about  10  eV 

energy  is  capable  of  ionizing  the  upper  states  of  any  nitrogen  ion.  Thus, 

this  relatively  weak  line  has  also  reduced  the  N3+  population  from  about 
20%  to  about  13%  and  increased  the  N54  population  from  about  9%  to  18%. 

In  Figure  19  a  line  at  46.556  eV  connecting  states  2  and  3  has 
been  added.  This  line  has  the  expected  effect  of  bringing  state  3  into 
approximate  equilibrium  with  states  1  and  2  but  more  importantly,  it  is 
capable  of  ionizing  any  state  of  N04  or  N14  and  all  but  the  lowest  states 
of  N24,  N3+,  and  N44.  The  steady  state  configuration  for  this  case  is: 

0+  1+  24  34  4+  5+  64 

7xl0-19  4x10" 11  2x10" 7  10“4  1.6x10“ 3  .998  8xl0“23 

The  effect  has  been  to  raise  the  mean  state  of  ionization  from 
about  4  to  5.  Under  borderline  conditions,  radiation  will  virtually  cease 
from  affected  regions.  Both  N3+  and  N4+  are  good  radiators,  whereas  N54 
is  an  exceptional ly  poor  one. 


1 ,, 
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Figure  20  shows  the  barely  noticeable  effect  of  adding  the 
56.554  eV  line  connecting  states  1  and  4  to  the  previous  3  lines.  All 
states  below  state  4  are  reduced  in  population  by  about  a  factor  of  2  and 
no  other  discernable  event  occurs. 


Finally,  Figure  21  shows  a  comparison  of  the  case  with  no  radia¬ 
tion  to  that  with  all  four  lines  over  the  full  range  of  eigenstates.  It 
shows  that  low  ionization  states  are  depopulated  by  up  to  12  orders  of 
magnitude  by  the  radiation  (either  the  2-3  line  or  the  1-4  line  would 
accomplish  most  of  this),  and  the  upper  states  of  N3*  are  brought  into 
equilibrium  with  the  lower  states  of  N4+. 

The  most  important  meaning  of  the  series  of  calculations  is  hid¬ 
den  in  the  millimeter  increase  in  the  ground  state  of  N ^  at  269  eV, 
bringing  it  from  9%  to  99%  of  the  population,  and  the  3  millimeter  de¬ 
creases  in  the  ground  states  of  N 3+  and  N4+,  reducing  them  from  major  ions 
to  obscurity.  Historically,  cases  similar  to  this  have  been  calculated 
assuming  line  radiation  drains  energy  from  the  system  but  does  not  appre¬ 
ciably  affect  state  of  ionization.  The  set  of  cases  calculated  here  do 
not  follow  the  dynamic  behavior  of  shocked  material  and  so  a  firm  conclu¬ 
sion  cannot  be  drawn,  but  certainly  they  call  into  question  the  previous 
method  of  calculation. 
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